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INTRODUCTION 


The  model  of  longwave  propagation  developed  at  the  Naval  Ocean  Systems  Center 
(NOSC)  is  based  on  a  waveguide  mode  formulation.  To  determine  signal  levels  in  this 
approach  the  basic  problem  is  to  obtain  the  modal  solutions  to  the  specific  waveguide 
under  consideration  (Pappert  et  al.,  1970.  Morfitt  and  Shellman,  1976).  where 
complicated  propagation  paths  are  divided  into  horizontally  homogeneous  segments. 
The  parameters  of  the  segments  are  determined  by  the  earth’s  ground  conductivity,  the 
magnitude  and  orientation  of  the  geomagnetic  field  with  respect  to  the  direction  of 
propagation,  and  the  state  of  the  ionosphere.  Well-known  computer  programs  which 
make  the  necessary  calculations  for  a  single  set  of  propagation  parameters  are 
"WAVEGUID"  (Pappert  et  al..  1970)  and  "MODESRCH"  (Morfitt  and  Shellman.  1976). 
The  WAVEGUID  computer  program  and  related  programs  are  described  in  a  series  of 
reports  and  familiarity  with  the  important  elements  of  this  series  is  assumed  (Pappert. 
Gossard  and  Rothmuller,  1967;  Sheddy  et  al.,  1968;  Pappert.  Moler  and  Shockey.  1970; 
Morfitt  and  Shellman,  1976).  This  report  describes  a  modified  version  of  that  program 
designed  to  supply  calculations  of  long  wavelength  propagation  along  segmented 
propagation  paths.  This  program  is  called  the  Segmented  Waveguide  (SW). 

The  computer  program  obtains  waveguide  mode  solutions  for  very  low  frequencies 
and  low  frequencies  (VLF/LF).  The  program  allows  for  multiple  homogeneous 
segments  to  be  specified,  allowing  for  consideration  of  variations  in  the  earth-ionosphere 
waveguide.  Path  geometry  and  geophysical  parameters  can  be  computed  by  the 
program.  Ionospheric  disturbances  due  to  man-made  or  naturally  occurring  events  can 
also  be  modeled  using  the  program. 

Essential  features  of  this  program  include: 

•  Automatic  segmentation  of  the  propagation  path 

•  Allowance  for  presegmentation  of  the  propagation  path 

•  Allowance  for  variation  of  the  ionosphere  along  that  path 

The  program  operates  on  propagation  paths  defined  by  a  transmitter  location  and 
either  a  direction  or  a  receiver  location.  A  propagation  path  is  defined  as  the  great 
circle  on  a  spherical  earth.  Variation  of  the  geophysical  parameters  are  to  be  expected 
along  realistic  paths.  The  diurnal  condition  in  large  part  determines  the  significance  of 
the  other  parameters.  For  instance,  under  daytime  conditions  the  effect  of  variation  of 
the  geomagnetic  field  along  a  path  is  usually  not  significant  to  the  resulting  mode 
parameters.  The  program  incorporates  routines  for  calculating  the  parameters  of  the 
geomagnetic  field  and  for  selecting  the  ground  conductivity  at  any  point  on  the  earth’s 
surface. 

Each  of  the  path  segments  is  treated  as  a  horizontally  homogeneous  planar 
waveguide.  Earth  curvature  is  introduced  by  use  of  a  modified  refractive  index.  A  set  of 
possible  solutions  to  the  waveguide  mode  equation  must  be  input.  Each  of  these 
solutions  is  processed  by  an  iteration  routine.  Each  iteration  requires  computation  of 
ionospheric  and  ground  reflection  coefficients.  Calculation  of  the  ionospheric  reflection 
coefficients  requires  integration  of  the  coefficients  through  the  ionosphere.  Ai 
approximate  formulation  may  be  used  which  requires  a  secondary  set  of  complex  angles 
be  specified  by  the  user.  In  that  case,  the  ionospheric  reflection  coefficients  are 
calculated  for  the  secondary  set  of  angles.  These  coefficients  are  used  to  interpolate 
the  ionospheric  reflection  coefficients  during  the  iteration  of  the  primary  set  of  possible 
solutions.  This  interpolation  procedure  requires  much  less  computation  time  than  does 
the  more  exact  procedure. 
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The  set  of  solutions  for  the  first  homogeneous  segment  must  be  input  by  the 
user.  The  program  uses  the  results  of  up  to  three  consecutive  segments  to  extrapolate 
the  solutions  for  successive  segments.  This  reduces  the  number  of  iterations  which  are 
required  for  subsequent  segments  and  allows  for  the  tracing  of  mode  solutions  through 
a  wide  range  of  path  variations. 

The  primary  output  of  "SW"  is  data  which  may  be  used  in  mode  summing 
programs.  The  strength  of  the  electromagnetic  field  along  the  path  can  be  obtained 
with  either  of  two  mode  conversion  models,  one  denoted  "FULLMC"  (Pappert  and 
Shockey,  1972)  and  the  other  denoted  "FASTMC"  (Ferguson  and  Snyder,  1980). 
"FULLMC"  does  integrations  through  the  ionosphere  prior  to  calculating  mode 
conversion  coefficients  and  can  be  quite  slow  in  execution  time  whereas  "FASTMC” 
avoids  the  transionospheric  integrations  by  use  of  approximations  and  runs  quite 
quickly. 

PROGRAM  CONTROL 

Program  execution  is  controlled  by  strings  containing  mnemonic  words.  These 
strings  and  the  sequence  of  operations  initiated  by  them  are  described  below.  These 
control  strings,  as  well  as  variable  names  and  names  of  subroutines,  appear  in  upper 
case.  For  clarity,  the  control  strings  are  enclosed  in  single  quotes  and  names  of 
subroutines  are  enclosed  in  double  quotes.  Table  1  summarizes  these  control  strings. 


Table  1.  Summary  of  control  strings. 

ID  run  identification 

NAME  general  NAMELIST  data 

EIGEN  alternate  EIGEN  list  input 

PROFILE  i  ionospheric  specie  profiles  (i  is  1  or  3) 

COLFREQ  ionospheric  collision  frequency  profiles 

COORD  automatic  segmentation  of  propagation  path 

PRESEG  presegmentated  propagation  path 

QUIT  end  of  job 


’ID'  indicates  that  the  next  line  of  data  to  be  read  is  a  general  identification  for 
the  path  under  consideration.  This  identification  appears  in  the  printout.  It  is  also 
written  to  the  mode  parameter  output  file. 

’NAME’  initiates  the  reading  of  general  program  data  via  the  NAMELIST: 
DATUM.  Table  2  lists  all  variables,  their  type  (I  indicates  integers,  R  indicates  real 
(floating  point)  variables,  C  indicates  complex  variables),  their  units  where  applicable, 
and  their  initial  values.  If  a  variable  defines  an  array,  then  the  dimension  of  the  array 
follows  the  name  in  parentheses.  The  NAMELIST  input  format  is  quite  flexible  but 
requires  that  column  1  be  blank.  Variable  names  are  followed  by  an  equal  sign  and 
then  by  the  value  of  the  variable.  Array  variable  names  may  be  followed  by  a  string 
of  values  separated  by  commas  and/or  spaces.  Embedded  blanks  are  not  allowed  in 
the  variable  names.  Variable  types  must  be  considered;  for  example,  values  for  integer 
variables  may  not  contain  decimal  points,  but  values  for  real  variables  do  not  have  to 
have  decimal  points.  Logical  variables  may  be  specified  with  any  of  ’.true.’  ’.t.’  ’t’ 
’.false.’  ’.f.’  ’f’.  The  values  of  character  variables  must  be  enclosed  within  quotes. 
The  first  record  of  the  NAMELIST  input  must  contain  ’  &name  ’  where  ’name'  is  the 


NAMELIST  name  (in  this  case,  DATUM).  The  end  of  the  NAMELIST  is  indicated  by 
’&end’.  Some  of  the  values  to  be  found  in  the  following  text  are  of  the  form  A(N)  to 
indicate  A  to  the  power  N. 

’EIGEN’  allows  for  input  of  the  list  of  trial  solutions  to  be  made  from  a  file.  The 
format  of  the  input  is  the  same  as  with  NAMELIST.  The  control  string  is  followed 
by  the  name  of  the  file  (starting  in  column  9)  containing  the  NAMELIST  data.  The 
source  of  these  input  solutions  could  be  a  previous  run  with  the  program  or  from  one 
of  the  automatic  mode  searching  programs  such  as  the  one  described  by  Morfitt  and 
Shellman  (1976). 

'PROFILE  i'  initiates  reading  of  the  ionospheric  charged  particle  profile  data  used 
to  model  the  upper  boundary  of  the  earth-ionosphere  waveguide.  The  value  of  i 
indicates  the  number  of  ionosperic  species  to  be  used  and  it  must  have  one  of  two 
possible  values:  1  is  for  electrons  only  and  3  is  for  electrons  and  ions.  If  i  is  not 
specified,  a  value  of  1  is  assumed.  The  use  of  i  is  shown  below.  The  ’PROFILE  i’ 
string  is  immediately  followed  by  a  single  line  of  identification  for  the  profile.  The 
profile  is  input  starting  at  the  top  of  the  ionosphere  using  a  formatted  input.  Each 
line  contains  the  height  in  km,  the  electron  density  at  that  height  in  electrons  per 
cubic  centimeter,  and  if  ions  are  to  be  considered,  the  positive  ion  density  at  that 
height  in  ions  per  cubic  centimeter.  The  height  is  in  columns  1-7  and  the  electron 
and  ion  densities  are  required  to  be  in  columns  14-21  and  24-31,  respectively.  The 
end  of  the  profile  is  indicated  by  a  dummy  height  with  value  less  than  zero.  A 
maximum  of  50  heights  may  be  used.  If  i  is  1,  then  only  the  electron  density  is  read. 
Consequently,  only  the  electron  density  need  be  present  in  the  data.  If  i  is  3,  then 
the  electron  and  positive  ion  densities  are  read  and  the  negative  ion  density  is 
computed  by  subtracting  the  electron  density  from  the  positive  ion  density  (to  preserve 
charge  neutrality). 

In  the  integration  of  the  reflection  elements  through  the  ionosphere  the  program 
interpolates  exponentially  between  input  heights.  The  profile  should  contain  sufficient 
data  to  define  the  ionospheric  structure  with  height.  For  example,  an  exponential 
profile  should  consist  of  only  the  top  and  bottom  heights  and  densities.  Many 
regularly  spaced  heights  tend  to  slow  the  integration. 

A  purely  exponential  conductivity  profile  (electrons  only)  may  be  input  via  the 
NAMELIST  variables  BETA  and  HPRIME. 

Additional  specie  parameters  are  needed  for  the  waveguide  mode  computations  and 
are  described  below. 

’COLFREQ’  initiates  reading  of  an  ionospheric  collision  frequency  profile.  This 
allows  use  of  nonexponential  collision  frequencies.  The  ’COLFREQ’  string  is 
immediately  followed  by  the  collision  frequency  profile,  starting  with  the  highest  height 
and  ending  with  a  dummy  height  of  value  less  than  zero  as  in  the  case  of  specie 
profile  described  above.  These  heights  need  not  match  those  used  under  ’PROFILE  i’. 
The  format  of  the  data  is  the  same  as  used  for  ’PROFILE  i’  except  rh.t  cc!ii«Jcr 
frequencies  for  all  species  must  be  input  since  the  negative  ion  col  isi  *'  f'er1ue~  v 
cannot  be  simply  computed  from  the  other  two.  If  only  electrons  are  being  used,  then 
only  that  collision  frequency  need  be  present.  As  with  the  specie  profiles,  the  program 
interpolates  exponentially  between  input  heights. 

An  exponential  collision  frequency  specification  may  be  input  via  the  NAMELIST 
variables  EXPNU  and  COEFNU. 


3 


'COORD'  initiates  automatic  segmentation  of  the  propagation  path.  This  string 
must  be  placed  after  all  pertinent  data  have  been  read.  This  option  is  best  applied  to 
simple  cases  suc'i  as  all  daytime.  The  basic  input  consists  of  the  path  specification, 
the  environment,  and  the  starting  mode  solutions. 

’PRESEG*  allows  for  previously  determined  segments  to  be  used  along  the 
propagation  path.  This  option  requires  most  of  the  same  inputs  as  ’COORD’  except 
that  the  user  supplies  the  distances  at  which  segments  begin.  At  each  segment  the 
user  has  the  option  of  specifying  the  parameters  of  the  geomagnetic  field,  the  ground 
and  the  ionosphere. 


CALCULATION  OF  MODE  PARAMETERS 

The  inputs  to  the  mode  equation  computations  are  supplied  by  geophysical 
routines  and/or  by  the  user.  The  subroutine  which  controls  the  calculations  is 
"WVGUID".  Most  user  supplied  data  values  are  input  to  the  program  via  NAMELIST. 
These  parameters  are  summarized  in  table  2  and  are  discussed  in  more  detail  below. 
In  table  2  the  data  types  are  Integer  (I),  Real  (R),  and  Complex  (C). 


Varible 

Type  Default 

FREQ 

R 

0.0 

RHO 

R 

0.0 

AZIM 

R 

0.0 

CODIP 

R 

0.0 

MAGFLD 

R 

0.0 

SIGMA 

R 

4.64 

EPSR 

R 

81.0 

BETA 

R 

0.0 

HPRIME 

R 

0.0 

TLONG.TLAT  R 

0.0.0.0 

RBEAR 

R 

720.0 

RLONG.RLAT  R 

O.O.O.O 

DRMIN 

R 

0.125 

DRMAX 

R 

0.5 

DM  AX 

R 

20.0 

YEAR 

1 

0 

MONTH 

1 

0 

DAY 

1 

0 

GMT 

R 

0.0 

NPRINT 

1 

1 

NPROF 

1 

1 

Table  2.  NAMELIST  Inputs. 

Description 
Frequency  in  kHz 

Distance  from  the  transmitter  of  the  current  segment 
in  Mm. 

Magnetic  azimuth  angle  in  degrees  east  of  magnetic 
north. 

Magnetic  co-dip  angle  in  degrees. 

Intensity  of  the  geomagnetic  field  in  Webers/square 
meter. 

Ground  conductivity  in  Siemens/meter. 

Dielectric  constant  of  the  ground. 

Slope  of  the  exponential  profile  in  km(-l) 

Reference  height  of  the  exponential  profile  in  km. 

Transmitter  coordinates  in  degrees 
Geographic  bearing  of  the  path  in  degrees 
Receiver  coordinates  in  degrees 

Minimum  distance  step  size  in  Mm 
Maximum  distance  step  size  in  Mm 
Maximum  distance  in  Mm 

Year 

Month  number  with  January  being  1 

Day  of  the  month 

Greenwich  meridian  time  in  hours 

Flag  used  to  control  the  amount  of  print  out 

Flag  to  indicate  which  form  of  the  profile  is  to  be 

used 


MIDPNT 


I  0 


Flag  to  indicate  that  only  the  midpoint  ic  to  tc 
considered 

IGCD  I  0  Flag  to  indicate  that  the  computed  distance  between 

the  transmitter  and  receiver  is  to  be  used 

IGND  I  0  Flag  to  indicate  that  the  ground  conductivity  is  to  be 

determined  from  the  ground  map 

MDIR  I  0  Flag  used  to  reverse  the  direction  of  the  magnetic  field 

EIGEN(60)  C  0.0, 0.0  Initial  solutions  for  the  waveguide  modes  in  degrees. 

TLIST(30)  C  0.0, 0.0  Angles  where  the  reflection  coefficients  are  computed 

for  the  inexact  interpolation  routine  in  degrees 
DTHETA  C  0.01,0.001  Change  in  the  mode  solution  used  to  define  the 

iteration  in  degrees 

LUB  C  0.05,0.005  Tolerance  used  to  test  the  differential  change  in  the 

mode  solution.  Used  to  stop  the  iteration 
OEIGEN  C  0.05,0.005  Tolerance  used  to  define  duplicate  modes  in  degrees 

FTOL  R  1.0  Tolerance  used  to  determine  if  the  mode  equation  has 

been  satisfied  by  the  solutions 

THTINC  R  0.5  Maximum  change  in  degrees  of  either  the  real  or 

imaginary  part  of  the  mode  solution  from  one  iteration 
to  the  next 

MAXITR  I  7  Maximum  number  of  iterations  to  attempt 

ALPHA  R  3.14xl0(-4)  The  earth  curvature  correction  factor  in  km(-l) 

H  R  50.0  Height  at  which  the  modified  refractive  index  is  unity 

in  km 

D  R  0.0  Height  at  which  the  integration  through  the 

ionosphere  is  stopped  in  km 

PREC  R  2.0  Factor  which  controls  the  precision  of  the  reflection 

coefficient  integration 

WR0  R  2.5x10(5)  Value  of  omega  sub  r  used  to  define  the  reference 

height 

ATNMAX  R  50.0  Maximum  attenuation  rate  of  modes  to  be  retained  in 

dB/Mm 

DEBUG  I  0  Flag  used  to  generate  additional  printout  for  debugging 

purposes 

TYPITR  I  0  Flag  used  to  define  the  form  of  the  mode  equation 

RPOLY  I  1  Flag  used  to  define  the  reflection  coefficient  calculation 

NRTLST  I  5  Number  of  points  to  use  in  the  interpolation  of  the 

reflection  coefficients  during  inexact  iterations 
LUNIT7  I  7  Logical  unit  number  to  which  the  mode  parameters 

data  are  output 

CHARGE(3)  R  -1.0,  Charge  of  the  ionospheric  species 

1.0.-1.0 

MRATIO(3)  R  1.0,  Ratio  of  the  mass  of  the  ionospheric  spcr'<=.  to 

2*5.8x10(4)  that  of  electrons 

COEFNU(3)  R  1.816x10(11)  Collision  frequency  of  the  ionospheric  species  at 
2*4.54x10(10)  the  ground  in  collisions/sec 
EXPNU(3)  R  -0.15,  Exponential  slope  of  the  collision  frequency 

2*-0.15  in  km(-l) 


The  radio  frequency  in  kHz  is  specified  by  the  variable  FREQ.  The  input  to 
subroutine  ’’WVGUID”  includes  an  ionospheric  profile.  The  variable  NPROF  controls 
which  ionospheric  profile  is  to  be  used.  A  value  of  0  indicates  that  the  profile  input 
via  ’PROFILE  i’  is  to  be  used.  A  value  of  1  indicates  that  an  exponential  electrons 
only  profile  is  to  be  used.  The  profile  is  specified  by  the  exponential  slope  BETA  in 
km(-l)  and  a  reference  height  HPRIME  in  km  (Wait  and  Spies,  1964).  A  value  of  2 
for  NPROF  is  used  to  indicate  that  a  series  of  profiles  will  be  input.  This  option  only 
applies  to  ’PRESEG’.  The  profiles  are  to  be  input  using  the  same  format  as  described 
under  ’PROFILE  i',  including  the  control  string.  There  must  be  a  profile  for  each 
segment. 

Additional  specie  parameters  are  needed.  The  number  and  order  of  these  specie 
parameters  must  be  consistent  with  the  charged  particle  densities  of  ’PROFILE  i'.  The 
charges  of  the  species  are  input  as  CHARGE  (i.e.,  CHARGE  =  -1,  +1,  -1).  The 
masses  of  the  species  relative  to  that  of  an  electron  are  input  as  MRATIO.  The 
collision  frequencies  may  be  defined  with  ’COLFREQ’  (nonexponential)  or  with  the 
variables  COEFNU  in  collisions  per  second  and  EXPNU  in  km(-l).  The  collision 
frequency  v  at  an  altitude  z  is  then  defined  by 

v  =  COEFNU  *  exp(EXPNU  *  z) 
where  z  is  in  km. 

Parameters  of  the  geomagnetic  field  are  specified  by  AZIM.  the  angle  between 
magnetic  north  and  the  direction  of  propagation  in  the  horizontal  plane  measured  in 
degrees  east  of  north,  CODIP,  the  magnetic  co-dip  angle  measured  from  the  vertical 
(i.e.,  the  north  pole  has  a  CODIP  of  0),  and  MAGFLD,  the  magnetic  intensity  in 
Webers  per  square  meter  or  in  Gauss.  The  magnitude  of  MAGFLD  is  tested.  If  it  is 
greater  than  10(-2)  then  the  input  value  is  assumed  to  be  in  Gauss  and  is  mu'^iplied 
by  10(-4). 

MDIR  is  a  flag  that  when  set  to  1,  causes  the  direction  of  propagation  as  input 
to  be  reversed.  This  allows  for  development  of  a  data  set  appropriate  to  examining 
transmitter  deployment. 

Ground  conditions  are  specified  by  SIGMA,  the  conductivity  in  Siemens/m,  and 
EPSR.  the  relative  dielectric  constant. 

The  presegmentation  option  allows  the  magnetic  field  and  ground  parameters  to  be 
varied  by  the  user. 

The  correction  for  earth  curvature  is  controlled  by  ALPHA  in  km(-l)  which  is 
defined  as  2  over  the  radius  of  the  earth.  For  a  curved  earth  ALPHA  is  3.14xl0(-4) 
km(-l)  and  for  a  flat  earth  ALPHA  is  0. 

Ionospheric  altitude  parameters  are  H,  which  is  the  height  in  km  at  which  the 
modified  refractive  index  is  unity  and  is  the  height  to  which  the  mode  solutions  are 
referenced;  D,  which  is  the  height  in  km  below  which  ionospheric  effects  can  be 
ignored.  D  must  be  equal  to  or  greater  than  the  bottom  height  of  the  ionospheric 
profiles  used.  It  is  usually  sufficient  to  choose  H  equal  to  D.  The  choice  of  D  and 
H  is  also  discussed  by  Pappert  et  al.  (1967). 


The  trial  eigenangles  follow  the  variable  name  EIGEN  which  is  a  complex  variable. 
Up  to  60  eigen  angles  may  be  input.  If  little  is  known  about  the  expected  solutions 
for  a  given  set  of  conditions,  a  set  of  approximate  solutions  may  be  obtained  using  a 
TL1ST.  The  TLIST  is  a  list  of  as  many  as  30  complex  angles  which  are  used  to  set 
up  an  interpolation  matrix  of  the  ionospheric  reflection  elements  (Sheddy  et  al..  1968i. 
The  program  then  uses  this  matrix  to  interpolate  reflection  elements  during  the  ite'atur 
process  used  to  obtain  mode  solutions.  These  solutions  are  referred  to  as  "inexacts" 
in  order  to  distinguish  them  from  the  more  accurate  solutions  using  integrated  reflection 
coefficients.  The  variable  NRTLST  determines  the  maximum  number  of  TLIST  angles 
used  in  each  interpolation.  During  the  inexact  iteration  process,  the  program  computes 
the  magnitude  of  the  complex  difference  between  the  current  value  of  the  solution  and 
each  of  the  TLIST  angles.  The  program  orders  the  TLIST  angles  from  th<  -.,r:  .T t  : 
the  largest  difference  and  selects  the  first  NRTLST  of  them  to  be  used  in  the 
interpolation.  This  improves  the  accuracy  of  the  interpolated  reflection  coefficients  and 
reduces  the  number  of  terms  used. 

If  more  than  30  EIGENs  are  input,  then  the  program  sorts  the  angles  according  to 
their  attenuation  rate  and  deletes  those  with  attenuation  rates  greater  than  a  user 
specified  maximum.  The  initial  value  of  this  maximum  is  ATNMAX.  If  there  are  still 
more  than  30  angles,  then  the  maximum  attenuation  rate  is  reduced  by  5  dB  Mm  and 
the  input  list  is  sorted  again.  This  process  is  repeated  until  there  are  less  than  30 
angles  in  the  list. 

If  the  number  of  EIGEN  or  TLIST  inputs  varies  from  one  NAMELIST  to  the  next, 
then  each  EIGEN  or  TLIST  list  should  be  terminated  with  a  zero.  If  RPGLY  is  not  0 
and  the  first  value  of  TLIST  is  0,  then  TLIST  is  set  equal  to  the  first  30  EIGENs. 

The  Newton-Raphson  iteration  process,  used  to  find  the  eigenangles  which  su  v 
the  modal  equation  is  described  by  Sheddy  et  al.  (1968).  Iteration  is  performed  fur 
each  input  EIGEN.  The  iteration  stops  when  the  maximum  number  of  iterations 
(MAXITR)  is  exceeded  or  when  the  change  in  the  real  and  imaginary  parts  of  tin- 
solution  is  calculated  to  be  less  than  the  real  and  imaginary  parts  of  LUB,  respectively. 

The  type  of  solution  obtained  is  determined,  in  part,  by  RPOLY.  which  can  have 
three  values:  0  for  exact  solutions  only,  2  for  inexact  (approximate)  solutions  only,  and 
1  for  inexacts  computed  and  used  as  inputs  to  obtain  exact  solutions.  The  use  of 
RPOLY  equal  to  1  is  described  more  fully  below. 

The  flag  TYPITR  is  used  to  obtain  vertically  polarized  modes  only  (TYPITR  equal 
to  1)  or  horizontally  polarized  modes  only  (TYPITR  equal  to  2).  It  is  physically 
meaningful  to  apply  this  option  only  for  nearly  isotropic  conditions,  no  magnetic  field 
(MAGFLD  set  to  0).  or  east  to  west  and  west  to  east  propagation  at  the  geomagnetic 
equator  (CODIP  is  90  and  AZIM  is  90  or  CODIP  is  90  and  AZIM  is  270). 

To  ensure  consistent  mode  sums  and  eliminate  redundant  solutions,  each  exact  and 
inexact  EIGEN  solution  is  tested  for  several  conditions.  The  first  is  that  the  imaginary 
part  of  the  solution  must  be  less  than  zero  in  order  to  have  attenuating  rr  der.  Yt  .■> 
second  is  that  the  magnitude  of  the  modal  equation  must  be  less  than  F  S  ut.  This 
parameter  is  tested  only  for  final  solutions  and  only  if  the  number  of  iterations  required 
to  obtain  the  solution  is  greater  than  or  equal  to  MAXITR.  It  is  generally  true  that  it 
the  iteration  stops  because  the  change  in  the  mode  solution  is  less  than  LUB,  then  the 
value  of  the  mode  equation  is  small.  There  are  instances  in  which  the  test  on  FTOL 
will  still  fail.  Consequently,  the  value  of  FTOL  is  set  very  high  in  order  to  allow  th« 
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program  to  continue  execution.  In  some  cases  the  user  may  want  to  modify  the 
default  value  in  order  to  perform  special  tests.  If  RPOLY  is  1,  the  inexacts  are 
treated  as  intermediate  results.  The  third  test  is  that  the  value  of  the  EIGEN  solution 
must  be  different  from  all  previous  solutions  by  an  amount  DEIGEN  which  is  input  as 
a  complex  number.  The  real  and  imaginary  parts  of  DEIGEN  are  the  tolerances  for 
the  real  and  imaginary  parts  of  the  EIGEN  solutions,  respectively.  If  one  of  the  above 
tests  results  in  a  mode  being  dropped  from  the  list  of  solutions,  then  the  program 
follows  the  procedures  outlined  below  under  the  discussion  of  mode  tracing. 

Subroutine  "WVGUID”  computes  and  prints  attenuation  rate  in  dB/Mm,  phase 
velocity  relative  to  the  speed  of  light,  the  magnitude,  and  phase  of  Wait's  excitation 
factor  (Wait,  1962)  at  the  ground  in  dB  and  radians. 

The  headings  for  the  number  of  iterations  to  go  from  the  input  angle  to  the  final 
solution,  the  final  solution,  the  magnitude  of  the  modal  equation,  and  the  magnitude  of 
the  polarization  mixing  ratio  are  printed  as  ITER,  EIGEN,  MAG  F,  and  MAG  P, 
respectively.  The  attenuation  rate,  phase  velocity  relative  to  the  speed  of  light, 
magnitude,  and  phase  of  Wait’s  excitation  factor,  and  the  final  solution  references  to 
the  ground  are  printed  under  the  headings  of  ATTEN,  V/C,  WAIT’S  EXC,  and 
THETA’,  respectively. 

The  parameters  YEAR,  MONTH,  DAY,  and  GMT  are  used  only  to  pass  the  values 
to  the  output  files.  These  parameters  are  useful  for  helping  to  identify  the  output 
data  to  programs  which  may  use  this  information. 

PATH  CALCULATIONS 

These  calculations  are  controlled  by  ’’GCPATH’’.  They  can  be  divided  into  three 
classes.  The  first  automatically  computes  geometry  and  geophysical  parameters  which 
are  obtainable  from  just  the  location  of  a  point  on  the  propagation  path  or  uses 
presegmented  distances  and  geophysical  parameters.  In  addition,  the  program 
extrapolates  the  EIGEN  list  and  TUST  so  as  to  trace  modes  along  a  path. 

PATH  GEOMETRY 

All  geometry  calculations  are  for  a  spherical  earth.  Inputs  to  this  portion  of  the 
program  consist  of  transmitter  and  receiver  locations,  path  length,  and  path  increments. 
Transmitter  and  receiver  longitude  and  latitude  are  input  with  TLONG  and  RLONG  and 
TLAT  and  RLAT.  The  convention  used  in  the  program  is  east  longitude  and  south 
latitude  are  negative.  An  alternate  input  for  the  receiver  position  is  its  bearing, 
RBEAR,  in  degrees  east  of  north.  In  the  execution  of  the  program  RBEAR  is  tested. 
If  RBEAR  is  720,  then  the  program  uses  RLONG  and  RLAT  to  define  the  path.  If 
RBEAR  is  not  720,  then  it  and  the  input  path  length,  DMAX,  are  used  to  define  the 
path.  DMAX  is  specified  in  megameters  and  must  be  less  than  or  equal  to  20.  If 
RLONG  and  RLAT  are  used,  there  are  two  path  lengths  possible.  If  the  parameter 
IGCD  is  equal  to  1,  the  path  length  is  set  equal  to  the  computed  short  great  circle 
distance  between  the  transmitter  and  receiver.  If  IGCD  is  equal  to  0,  the  path  length 
is  unchanged  from  what  was  input  (DMAX).  If  the  path  bearing  is  input,  the  path 
length  is  always  DMAX. 

The  starting  value  of  the  distance  from  the  transmitter  is  input  as  RHO  and  is  in 
Mm.  The  path  increments  are  controlled  by  the  mode  tracing  results  and  the  variables 
DRMIN  and  DRMAX  which  are  all  in  units  of  Mm.  The  procedures  of  the  values  are 
described  as  follows. 
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The  geomagnetic  field  is  computed  at  the  first  path  point  defined  by  K.i u  r-nc.  _> 
the  beginning  of  each  path  segment.  The  ground  conductivity  and  relative  dielectric 
constant  are  specified  by  the  user  through  SIGMA  and  EPSR  via  NAMELIST  or  by 
searching  the  ground  map.  If  IGND  is  0.  the  ground  map  is  not  searched  and  ground 
conditions  are  assumed  constant,  as  input,  for  the  entire  path  for  ’COORD’  or  as 
varied  by  the  user  for  'PRESEG’.  If  IGND  is  1,  the  DECO-NRL  10  level  ground 
conductivity  map  (Hauser,  Garner,  and  Rhoads,  1969)  is  searched  for  the  appropriate 

values  of  SIGMA  and  EPSR  at  the  beginning  of  each  path  segment. 

If  the  value  is  assumed  that  the  entire  propagation  path  can  be  described 

adequately  by  the  conditions  at  the  midpoint  of  the  path,  then  the  path  conditions  at 

that  point  can  be  obtained  if  MIDPNT  is  set  to  1.  The  subsequent  modal  calculations 
will  then  be  for  the  midpoint  conditions. 

At  the  transmitter  and  at  the  end  of  each  path  segment,  the  parameters  to  be 
used  in  the  "WVGUID’’  calculations  are  printed  next  to  the  heading,  PROPAGATION 
PATH  PARAMETERS,  as  described  below.  In  addition,  the  distance  in  megameters 
from  the  transmitter,  the  coordinates  and  the  geographic  bearing  of  the  path  at  the 
current  point  are  printed  under  the  headings  RHO,  LAT,  and  BEAR,  respectively.  If 
the  midpoint  option  is  being  used,  then  the  above  information  for  the  midpoint  is 
printed. 

PRESEGMENTATION 

In  some  instances,  user  segmentation  of  the  propagation  path  is  desired.  The 
control  string  ’PRESEG’  allows  arbitrary  segmentation  of  the  path.  This  is 
accomplished  by  a  succession  of  data  lines  in  list  directed  format  containing  values  for 
path  distance  in  Mm.  AZIM.  CODIP.  MAGFLD,  SIGMA.  EPSR.  BETA,  and  HPRIME, 
respectively.  List  directed  input  is  accomplished  by  entering  values  separated  by 
commas  or  spaces.  There  must  a  data  entry  for  each  request  in  the  input  list.  If  a 
value  is  not  to  change  from  one  data  line  to  the  next,  then  the  value  need  not  be 
entered  but  its  omission  must  be  indicated  by  a  pair  of  commas.  The  first  value  of 
path  distance  need  not  be  zero.  The  presegmentation  is  terminated  by  a  distance 
value  of  40. 

If  NPROF  is  0,  the  ionospheric  profile  is  constant  for  the  path  and  is  defined  by 
’PROFILE  i’.  If  NPROF  is  1,  then  the  values  of  BETA  and  HPRIME  on  the 
presegmentation  data  lines  are  examined.  If  BETA  is  zero,  then  the  previously  defined 
values  of  BETA  and  HPRIME  are  used.  The  latter  may  be  input  via  NAMELIST  so 
that  a  constant  ionosphere  for  the  path  can  be  obtained  by  using  NAMELIST  input. 
If  either  of  the  values  of  BETA  and  HPRIME  are  to  change,  both  must  be  input.  If 
NPROF  is  2,  then  each  presegmentation  data  line  must  have  a  corresponding  ’PROFILE 
i’  profile  specification  on  logical  unit  3. 

The  following  conventions  are  used  for  using  the  values  of  the  presegmentation 
data.  If  the  value  of  MAGFLD  is  zero  or  blank,  then  the  magnetic  field  parameters 
are  calculated.  If  a  nonzero  value  is  specified,  then  all  of  the  magnetic  field 
parameters  are  taken  from  the  presegmentation  data.  If  constant  magnetic  p.j.  r,>  tr;  ■_ 
are  desired  along  the  path,  the  values  must  be  specified  on  each  presegmentation  data 
line.  If  the  value  of  SIGMA  is  not  entered,  then  the  previously  specified  value  of 
SIGMA  and  EPSR  are  used.  Otherwise,  the  values  of  these  parameters  are  taken  from 
the  presegmentation  data  line.  Constant  values  for  the  whole  path  may  be  specified 
via  NAMELIST  or  in  the  first  presegmentation  card.  If  IGND  is  1,  then  the  ground 


map  is  searched  and  the  values  of  SIGMA  and  EPSR  on  the  cards  are  ignored,  if 
BETA  is  not  entered,  the  currently  defined  values  are  used  for  the  electron  density 
profile.  Even  if  only  one  value  in  the  pairs  SIGMA,  EPSR,  and  BETA,  HPRIME  is  to 
be  changed,  both  values  must  be  specified. 

MODE  TRACING 

Efficient  computation  of  mode  parameters  along  the  propagation  path  is  best 
achieved  by  using  RPOLY  set  to  1,  which  will  be  assumed  for  the  rest  of  this 
discussion.  At  the  first  point  on  the  path,  solutions  are  best  obtained  by  using  a 
TLIST  composed  of  angles  which  are  believed  to  be  approximately  correct  and  an 
EIGEN  list  of  many  regularly  spaced  angles  such  as  88,  -1.  87,  -1,  86,  -1,  etc. 
Alternatively,  the  EIGEN  list  should  be  the  list  of  approximately  correct  solutions  with 
the  TLIST  set  to  zero  or  the  ’EIGEN'  control  string  could  be  used  to  specify  solutions 
from  some  other  source.  The  program  computes  inexact  solutions  for  the  conditions  at 
the  first  point  on  the  path.  After  exhausting  the  EIGEN  list,  obtaining  inexact 
solutions,  and  deleting  of  those  failing  the  tests  discussed  above,  exact  solutions  are 
computed  using  the  results  of  the  inexacts. 

Now  the  discussion  must  be  separated  for  the  two-path  segmentation  options.  For 
the  ’COORD’  option,  the  second  point  on  the  path  is  DRMIN  from  the  transmitter. 
For  this  point,  the  final  solutions  for  the  first  point  are  placed  in  both  TLIST  and 
EIGEN  and  the  same  process  of  calculation  of  inexacts  and  exacts  is  repeated.  If 
DRMIN  is  not  too  far  from  the  transmitter  and/or  the  geophysical  parameters  do  not 
change  too  much,  then  this  step  in  the  extrapolation  process  is  quite  efficient.  Now 
the  program  has  two  sets  of  final  solutions  and  makes  a  linear  extrapolation  for  TLIST 
and  EIGEN  angles  for  the  third  point  on  the  path  which  is  twice  DRMIN  from  the 
transmitter.  The  sequence  of  calculations  for  the  inexact  and  exact  solutions  is 
repeated.  For  the  fourth  and  all  subsequent  points  on  the  propagation  path,  the 
program  uses  the  previous  sets  of  final  solutions  to  make  second  order  extrapolations 
for  TLIST  and  EIGEN  angles.  The  distance  increments  are  chosen  as  described  below. 

As  the  program  steps  out  along  the  propagation  path,  modal  solutions  may  be  lost 
or  removed.  First,  a  mode  may  be  lost  in  the  screening  process  in  subroutine 
"WVGUID"  as  described  above.  At  the  first  point  on  the  path,  modes  may  be 
overlooked  simply  because  of  lack  of  adequate  trial  solutions  and/or  more  than  one 
EIGEN  input  resulting  in  the  same  final  solution,  perhaps  due  to  closely  spaced  input 
EIGEN  values.  If  computations  are  being  made  at  the  transmitter  or  at  the  midpoint, 
it  is  acceptable  to  lose  a  solution  from  the  input  EIGEN  list.  At  all  other  points, 
when  a  mode  is  lost  execution  terminates  in  "WVGUID’’.  After  the  tests  on  the 
solutions  in  "WVGUID"  are  completed  at  the  first  point  on  the  path,  the  program 
assumes  that  it  has  a  complete  set  of  modes.  After  this  set  is  established,  solutions 
may  be  acceptably  removed  only  in  the  extrapolation  subroutine,  "EXTRAP".  The 
solutions  produced  by  "WVGUID”  for  the  current  segment  are  used  to  compute 
attenuation  rates.  Those  solutions  whose  attenuation  rate  exceeds  ATNMAX  are 
deleted  from  the  list.  The  location  of  the  solution  in  the  set  is  marked  and  its 
removal  is  indicated  by  a  blank  line  in  the  printout  of  solutions  produced  by 
subsequent  "WVGUID"  calculations. 

If  a  mode  is  lost  during  "WVGUID”  calculations,  the  path  point  is  moved  back  to 
about  halfway  between  the  current  point  (where  a  solution  was  lost)  and  the  previous 
point  (where  all  solutions  were  obtained).  The  actual  distance  depends  on  the  current 
value  of  the  distance  increment.  If  the  increment  is  greater  than  DRMIN,  then  the 


new  increment  is  chosen  that  it  is  an  integral  multiple  of  DRMIN  and  is  less  tha.'.  cr 

equal  to  half  the  previous  increment.  Geophysical  parameters  at  the  new  path  point 

are  computed,  the  EIGEN  list  is  revised  by  "EXTRAP",  and  "WVGUID"  calculations 
are  repeated. 

If  no  modes  are  lost  and  the  number  of  iterations  required  to  obtain  the  solutions 
is  less  than  or  equal  to  half  of  MAXITR,  then  the  distance  increment  is  increased  by 
DRMIN.  This  increase  in  the  distance  increment  can  continue  until  the  path  increment 
is  equal  to  DRMAX.  If  no  modes  are  lost  and  the  number  of  iterations  required  to 
obtain  the  solutions  is  greater  than  half  of  MAXITR,  the  path  increment  is  decreased 
by  DRMIN. 

If  modes  are  lost  and  the  separation  between  the  previous  point  (for  which  ail 
modes  were  found)  and  the  current  point  (for  which  modes  were  lost)  is  less  than  or 
equal  to  DRMIN,  then  the  distance  increment  is  halved.  The  geophysical  parameters 

for  the  new  path  point  are  linearly  interpolated  using  the  parameters  of  the  two  points 

at  which  geophysical  parameters  were  computed,  the  EIGEN  list  is  revised  by 
"EXTRAP”,  and  "WVGUID"  calculations  are  repeated.  Solutions  obtained  for 
interpolated  path  points  are  not  saved.  They  are  used  only  to  trace  the  mode 
solutions  between  the  points  for  which  the  geophysical  parameters  are  computed. 

If  no  modes  are  lost  and  the  number  of  iterations  required  to  obtain  the  solutions 
is  less  than  or  equal  to  half  of  MAXITR,  then  the  distance  increment  is  doubled. 
This  increase  in  the  distance  increment  can  continue  until  the  path  increment  is  equal 
to  the  distance  to  the  end  of  the  interpolation  interval.  If  no  modes  are  lost  and  the 
number  of  iterations  required  to  obtain  the  solutions  is  greater  than  half  of  MAXITR, 
the  path  increment  is  halved.  If  the  new  distance  increment  is  less  than  15  km,  then 
the  program  aborts. 

For  the  ’PRESEG’  option,  the  distance  increment  is  controlled  by  the  intervals 
between  the  presegmented  distances.  When  modes  are  lost,  the  interpolation  procedure 
for  cases  in  which  the  backup  interval  is  less  than  DRMIN  described  above  is  followed. 

OUTPUT 

Mode  parameters  from  the  program  are  written  to  the  logical  unit  whose  numerical 
value  is  LUNIT7.  The  first  line  of  data  ritten  contains  the  transmitter  location,  path 
bearing,  and  the  date  and  time,  as  input  through  NAMELIST.  The  identification  which 
followed  ’ID’  is  written  next.  If  no  identification  was  specified  with  ’ID’,  this  line  of 
data  is  blank.  The  identification  is  followed  by  a  sequence  of  lines  at  each  output 
distance. 

The  first  line  of  data  at  each  such  distance  contains  the  distance,  frequency, 
AZIM,  CODIP,  MAGFLD,  SIGMA,  EPSR,  and  the  reference  height  of  the  ionospheric 
profile.  In  descriptions  of  other  programs,  this  first  line  of  data  at  each  distance  is 
referred  to  as  the  RFACMSET  header.  This  header  line  is  followed  by  pairs  of  data 
lines,  one  pair  for  each  mode.  The  quantities  in  these  data  lines  •>'-  tue  rr .  !c 
solution  as  a  complex  angle  in  degrees,  a  flag,  Tl,  T2,  T3,  and  T1.  T'  0  o 

parameters  Tl,  T2,  T3,  and  T4  are  described  in  detail  by  Ferguson  and  Snyder  (1980). 
The  last  line  of  data  at  each  output  distance  is  blank.  These  data  ar*  suitable  for 
use  in  "FASTMC"  (Ferguson  and  Snyder,  1980). 
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If  the  program  fails  because  of  some  problem  at  the  first  path  point,  it  writes 
’Failure  at  RHO  1’  to  logical  unit  90.  Otherwise,  it  writes  the  distance  of  the  last 
point  for  which  "WVGUID"  successfully  completed.  If  the  end  of  the  path  is  reached, 
then  this  distance  is  output  as  40. 

SAMPLE  INPUT 

Sample  input  files  are  shown  in  figures  1  and  2.  In  the  first  sample  (figure  1), 
the  path  is  to  be  run  for  all  nighttime  conditions  assuming  all  seawater  ground.  The 
EIGEN  list  for  the  transmitter  is  input  directly  and  the  automatic  path  segmentation  is 
to  be  used. 

id 

Sample  run 
name 

icdatum  freq=23.4  h=50  d=75  Iunit7=7  atnmax=50 

tlong=150  tlat=20  rbear=10  dmax=10 

lub=.005  .0005  dtheta=.01  .001 

deigen=.05  .005  thtinc=.05 

beta=.43  hprime=87 

eigen  =  85.678  -0.206  84.595  -0.688  81.806  -0.609  81.027  -0.255 

77.653  -0.791  77.023  -0.269  73.199  -0.825  72.980  -0.300 

68.926  -0.264  68.599  -0.955  64.751  -0.213  63.907  -1.144 

60.457  -0.183 

&end 

coord 


Figure  1.  Sample  input  using  COORD  option. 

The  second  example  (figure  2)  is  a  much  more  complicated  case.  It  is  for  the 
same  path  of  the  first  sample,  but  the  ionosphere  is  to  be  varied  according  to  the 
diurnal  conditions  along  the  path  for  July  15  at  1612Z.  The  transition  from  night  to 
day  has  been  modeled  as  five  steps  starting  with  BETA  at  0.30  and  HPRIME  at  74 
and  ending  with  BETA  at  0.43  and  HPRIME  at  87.  In  addition,  the  ground 
conductivity  for  the  last  profile  takes  on  three  values:  4,  10(-2),  and  10(-3).  In  order 
to  improve  the  efficiency  of  the  mode  tracing,  additional  segmentation  has  been 
performed  so  that  each  ground  conductivity  is  processed  separately.  The  segmentation 
does  not  produce  final  output  that  is  monotonically  increasing  in  distance  from  the 
transmitter.  The  necessary  ordering  of  the  segments  must  be  performed  by  editing  the 
final  output  file  or  by  a  user  supplied  program.  The  initial  mode  solutions  for  each 
segment  have  been  already  calculated  and  stored  in  a  set  of  files  named 
XMTR202.MFx  where  x  ranges  from  0  to  6. 
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id 

Sample  run 
name 

&datum  freq=23.400  h=50.  Iunit7=7  atnmax=50. 
Iub=0.005  0.0005  dtheta=0.010  0.0010 
deigen=0.050  0.0050  thtinc=0.05 
year=84  month=  7  day=15  gmt=16.2 
tlong=  158.150  tlat=  21.417  dmax=4 
rbear=202.0  &end 
eigen  xmtr202.mf0 
preseg 

0.000,190.6,  50.8.0.350.4.E+00, 81., 0.30,74.0, 
40,0,0,0,0.0,0,0. 
eigen  xmtr202.mfl 
preseg 

0.500,190.8,  56.8, 0.337.4.E+00.81., 0.32, 76.2, 
40,0.0,0.0.0.0.0, 
eigen  xmtr202.mf2 
preseg 

0.960.190.9,  63.1.0.329.4.E+00.81..0.34.78.3. 
40.0.0.0.0.0.0,0. 
eigen  xmtr202.mf3 
preseg 

1.040,190.9,  64.2,0.327,4.E+00,81., 0.37.80.5, 

1 .240.190.9,  67.3.0.325.4.E+00.81., 0.37,80.5, 
40,0.0,0.0.0.0,0. 

eigen  xmtr202.mf4 
preseg 

1.340.190.9,  68. 8,0.324,4. E-i- 00,81. ,0.39,82. 7, 

1.540.190.9,  72.1,0.323.4. E+00, 81. ,0.39,82. 7, 
40,0,0,0.0.0.0.0. 

eigen  xmtr202.mf5 
preseg 

1.640.190.9,  73.8, 0.322, 4.E+00.81..0.41, 84.8, 

1.820.190.9,  76.8, 0.322, 4.E+00, 81. ,0.41, 84. 8, 
40.0,0.0,0,0.0,0, 

eigen  xmtr202.mf6 
preseg 

1.940.190.8,  78.9.0.322.4.E+00, 81., 0.43,87.0, 

2.120.190.8,  82.2,0.323,4.E+00,81. ,0.43, 87.0, 
2.300,190.7,  85.5,0.324,4.E+00,81. ,0.43, 87.0, 

2.480.190.6,  88. 8,0. 326, 4. E+00, 81. ,0.43, 87.0, 

2.660.190.6,  92.1.0.329.4.E+00, 81., 0.43,87.0. 
2.840,190.5,  95. 5,0.333,4. E+00, 81. ,0.43, 87.0, 
3.020,190.4,  98.9.0.337.4.E+00, 81. ,0.43, 87.0, 

40,0,0,0,0,0,0,0, 


Figure  2.  Sample  input  using  PRESEG  option. 
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SW:  SEGMENTED  WAVEGUID 


0001  c 
0002  c 
0003 

0004  1  c 

0005  1 

0006  1 
0007  1 

0008  1 
0009  1 

0010  1 
0011  1 
0012  1  c 

0013  1 

0014  1 

0015  1 

0016  1 
0017  1 

0018  1  c 

0019 

0020  1  c 

0021  1  c 

0022  1 
0023  1 

0024  1 

0025  1 

0026  1 
0027  1 

0028  1 
0029  1  c 

0030 

0031  1  c 

0032  1 

0033  1 

0034  1 

0035  1 

0036  1  c 

0037  1 

0038  1 

0039  1 

0040  1 

0041  1  c 

0042  1 

0043  1  c 

0044 

0045  1  c 

0046  1  c 

0047  1 

0048  1 

0049  1 

0050  1 

0051  1  c 

0052  c 
0053 
0054 
0055 
0056 
0057 


include  ’commonl .for/I ist’ 

common/ i nput/f  req , r ho , a  z i m , cod i p , magf I d , s i gma , epsr , beta , hpr i me , 

$  hprout 

common/path/pa thid,t long, tlat,r long, r I  at, rbear, dmax.drmi n, drmax, 

$  year, month, day ,gmt, npr i nt, nprof , npath, igcd, ignd,mdir, lost, 

S  I un i t7,  I x 

common/ i onosp/htl i st(50) , I n I i st(50, 3) , he  I i st(50) ,cf I i st(50,3) , 

S  charge(3) ,mratio(3) ,nrspec, Ihtmx, Ihtmn, lht,mbtmx,mhtmn,mht 

character*80  pathid 
integer  year, day 

rea I *4  f req, rho, az im, cod i p,magf I d, si gma, epsr, beta, hpr i me, hprout, 

S  1 1 ong, tl at, r long, r I  at, rbear ,dmax,drmi n,drmax,gmt, 

I  htl ist, Ini ist, he  I ist,cf I ist, charge, mratio 

include  ’commonl . i ni / I i st * 

initialize  commonl 

data  f  r eq/0 . / , rho/0 . / , az i m/0 . / , cod ip/0./, magf I  d/0 . / , 

S  sigma/4. 64/, epsr/81 ./, beta/0. /, hpr ime/0./, 

S  t I ong/0 . / , t I at/0 . / , r I ong/0 . / , r I at/0 . / , rbear/720 . / , 

S  dmax/20 . / , dr min/. 125/ , drmax/ . 5/ , md i r/0/, 

S  yea  r/0/ , month/0/ , day/0/ , gmt/0 . / , npr i nt/1 / , nprof /I / , 

S  igcd/0/, ignd/0/,mdir/0/, lunit7/7/, 

$  charge/-l . , 1 . , -1 ,/,mratio/l . , 2*58000. /,nrspec/l/ 

include  ’common2.for/l ist’ 

common/wg  in/el ist (2, 30) ,tl i st(2,30) , dtheta (2) , I ub(2) ,deigen(2) , 
S  th t i nc , f to  I , max i tr , a  I pha , h , d , prec , wrO , atnmax , debug , typ i tr , 

S  rpoly,nrtlst 

common/wg  out/tp(30) ,tterm(4,30) ,nterm(30) ,mode(30) , modes, nmds 

complex*8  tp, tterm.dthta 
integer  debug, typi tr, rpo ly 

rea 1*4  el i st, tl i st, dtheta , I ub,deigen,thti nc,ftol , a lpha,h,d,prec, 

S  wrO, atnmax 

equivalence  (dtheta, dthta) 

include  ’common2 . i n i / I i st ’ 

initialize  common2 

data  el ist/60*0./,tl ist/60*0./, 

$  dtheta/. 01, .001/, lub/.05, .005/,deigen/.05, .005/,thti nc/ .5/, 

$  f  to  I /1000 . / , max i tr/7/ , a  I pha/3 . 14e-4/ , h/50 . / , d/0 . / , prec/2 . / , 

S  wr0/2 . 5e5/ , atnmax/50 . / , debug , typ i tr/0, 0/ , rpo I y/1/ , nrt I  st/5/ 

name  I i st/da tum/f req, rho, az im,cod i p,magf I d, si gma , epsr , beta , hpr ime, 
S  t long, 1 1  at, r long, r I  at, rbear ,dmax,drmi n, drmax, 

S  year .month, day ,gmt, npr i nt, nprof ,mi dpnt, iged, i gnd ,md i r , 

S  I  uni t7, charge .mratio, coefnu.expnu, 

S  ei gen , 1 1 i st, dtheta , I ub.dei gen, thti nc,f to  I ,max i tr , 
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0058  $  a  I pha , h,d, prec, wrO, atnmax, typi tr ,  rpo ! y , nrt I st 

0059  c 

0060  complex  theta 

0061  character*  8  branch 

0062  character*40  fname 

0063  character*80  bed 

0064  c 

0065  dimension  coef nu(3) ,expnu(3) ,eigen(2,60) 

0066  c 

0067  c  Initialize  MAIN 

0068  data  nuf lag/0/, coefnu/1 .816ell, 2*4. 540e09/,expnu/3*- . 15/, 

0069  $  eigen/120*0./,midpnt/0/ 

0070  c 

0071  c  Unit:  Useage: 

0072  c  2  input  of  alternate  eigen  list 

0073  c  3  input  of  profiles  along  the  path 

0074  c  Iunit7  output  of  mode  parameters  along  the  path 

0075  c 

0076  c 

0077  10  read (5,1000, end=999)  bed 

0078  print  1001, bed 

0079  branch=bcd (1 : 8) 


0080 

if(branch  .eq. 

’id 

*  .or.  branch  .eq. 

’ID 

’) 

90 

to 

20 

0081 

if  (branch  .eq. 

’name 

’  .or.  branch  .eq. 

’NAME 

’) 

90 

to 

100 

0082 

if(branch  .eq. 

’eigen 

’  .or.  branch  .eq. 

’EIGEN 

’) 

go 

to 

130 

0083 

if  (branch  .eq. 

’prof i le 

’  .or.  branch  .eq. 

’PROFILE 

’) 

go 

to 

200 

0084 

if (branch  .eq. 

’coif req 

'  .or.  branch  .eq. 

’COLFREQ 

') 

go 

to 

250 

0085 

if  (branch  .eq. 

’preseg 

’  .or.  branch  .eq. 

’PRESEG 

’) 

go 

to 

400 

0086 

if  (branch  .eq. 

'coord 

’  .or.  branch  .eq. 

’COORD 

’) 

go 

to 

500 

0087 

if  (branch  .eq. 

’quit 

’  .or.  branch  .eq. 

’QUIT 

') 

go 

to 

999 

0088  c 

0089  print  *, ’ABORT  MAIN:  Control  card  not  recognized  ’ 

0090  stop 

0091  c 

0092  c  Path  identification 

0093  20  read (5, 1000)  path  id 

0094  print  1001, path  id 

0095  go  to  10 

0096  c 

0097  c  NAMELIST  input 

0098  100  read (5, datum) 

0099  if(nprint  .gt.  1)  print  datum 

0100  if(freq  .eq.  0.)  then 

0101  print  *, ’ABORT  MAIN:  FREQ  not  input  ’ 

0102  stop 

0103  end  if 

0104  if(magfld  .gt.  l.e-02)  magf ld=magf ld*l .e-04 

0105  go  to  10 

0106  c 

0107  c  Separate  EIGEN  list  input 

0108  130  read  (bed , 1004)  fname 

0109  open(unit=2,f i le=fname,status=’old’) 

0110  read(2, datum) 

0111  close(unit=2) 

0112  if(nprint  .gt.  1)  then 

0113  do  131  m=l,60 

0114  if (eigen (l,m)  .eq.  0.)  go  to  132 
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0115 

131 

e 

n 

E 

0116 

132 

print  1040,d, h, (e i gen (1 , k) , e i gen (2, k) , k=l , km) 

0117 

end  i  f 

0118 

go  to  10 

0119 

c 

0120 

c 

Prof i 1 e  i nput 

0121 

200 

read (bed, 1002)  number 

0122 

nrspec=maxO(l , number) 

0123 

nprof=0 

0124 

ca 1 1  prof i n(5,l ,50,npr i nt, nrspec, 1 htmx,htl ist, Ini ist) 

0125 

i f ( 1 htmx  . le .  0)  then 

0126 

print  *, ’ABORT  MAIN:  Ionospheric  profile  missing’ 

0127 

stop 

0128 

end  i  f 

* 

0129 

go  to  10 

0130 

c 

0131 

c 

Collision  frequency  profile  input 

0132 

250 

nuf 1 ag=l 

0133 

ca 1 1  prof i n (5,2, 50, npr i nt, nrspec,mhtmx, he  1 i st,cf 1 i st) 

0134 

if(mhtmx  . le.  0)  then 

0135 

print  *, ’ABORT  MAIN:  Collision  frequency  profile  missing’ 

0136 

stop 

0137 

end  i  f 

0138 

go  to  10 

0139 

c 

0140 

c 

Presegmented  path 

0141 

400 

npath=2 

0142 

go  to  600 

0143 

c 

0144 

c 

Automatic  path  segmentation 

0145 

500 

npath=mi dpnt 

0146 

c 

0147 

c 

Test  all  inputs  before  execution. 

0148 

c 

0149 

c 

Count  the  modes 

0150 

600 

do  602  m=l , 60 

0151 

if (eigen (l,m)  .eq.  0.)  go  to  603 

0152 

602 

nmds=m 

0153 

603 

i f (nmds  . 1 e .  0)  then 

0154 

print  *, ’ABORT  MAIN:  No  EIGEN  list  ’ 

0155 

stop 

0156 

end  i  f 

0157 

c 

Delete  duplicate  modes  using  DEIGEN 

0158 

if (nmds  .gt.  1)  then 

- 

0159 

m=l 

0160 

1=2 

0161 

605 

i  f (abs (ei gen (1 ,m) -ei gen (1 , 1 ) )  .It.  deigen(l)  .and. 

0162 

$  abs (e i gen (2,m) -e i gen(2, 1 ) )  .It.  deigen(2))  then 

0163 

c 

Found  a  match  so  drop  this  mode. 

0164 

do  607  k=l,nmds 

0165 

eigen(l,k)=eigen(l,k+l) 

0166 

607 

eigen(2,k)=eigen(2,k+l) 

0167 

e i gen  (1, nmds) =0. 

0168 

eigen(2,nmds)=0. 

0169 

nmds=nmds-l 

0170 

i f  ( 1  .  1 e .  nmds)  go  to  605 

0171 

end  i  f 
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% 

0172 

i f  ( 1  .It.  nmds)  then 

■’  t 

0173 

1  =  1+1 

{, 

0174 

go  to  605 

0175 

end  i  f 

0176 

if(m  .It.  nmds)  then 

0177 

m=m+l 

1 

0178 

1  =m+l 

0179 

go  to  605 

0180 

end  i  f 

0181 

end  i  f 

0182 

c 

0183 

610 

if (nmds  .gt.  30)  then 

0184 

c 

Too  many  modes  input,  reduce  the  number  by  deleting  input 

0185 

c 

eigen  list  values  which  have  attenuation  rates  in  excess 

0186 

c 

of  atnmax  and  re-count  the  modes 

0187 

capk=l . /(I .- . 5*alpha*h) 

0188 

aconst=-182 0426*f req 

0189 

atnmx=atnmax 

'*• 

0190 

611 

nm=0 

0191 

do  614  m=l,nmds 

0192 

if (eigen(l,m)  .eq.  0.)  go  to  615 

0193 

theta=cmp 1 x (e i gen (1 , m) , e i gen (2 , m) ) * ( . 01745329252 , 0 . ) 

0194 

i f (aconst*a imag (capk*cs i n (theta))  . le.  atnmx)  then 

0195 

if(nm  .eq.  30)  then 

4 

*>’■ 

0196 

antmx=atnmx-5 . 

0197 

go  to  611 

i.' 

0198 

else 

0199 

nm=nm+l 

0200 

el ist(l,nm)=eigen(l,m) 

0201 

el ist(2,nm)=eigen(2,m) 

, 

0202 

end  i  f 

- 

0203 

end  i  f 

0204 

614 

conti nue 

0205 

615 

nmds=nm 

* 

0206 

if(nprint  .gt.  1) 

0207 

$  print  1042,atnmax,  (el ist(l,k) ,el ist(2,k) ,k=l,nmds) 

*■ 

0208 

e  1  se 

* 

0209 

c 

Keep  all  input  modes. 

i 

0210 

do  616  m=l,nmds 

t' 

0211 

el  i st (1 ,m) =e i gen (1 ,m) 

0212 

616 

el ist(2,m)=eigen(2,m) 

0213 

end  i  f 

> 

0214 

if (nmds  .It.  30)  then 

4, 

0215 

el ist(l,nmds+l)=0. 

£ 

0216 

el i st(2, nmds+l)=0 . 

% 

0217 

end  i  f 

( 

0218 

c 

— 

0219 

i  f  ( r  po 1 y  .eq.  1  and.  1 1 i st (1,1)  .eq.  0.)  then 

0220 

do  619  m=l,nmds 

*• 

0221 

tl i st(l ,m)=el i st(l ,m) 

J« 

0222 

619 

tl ist(2,m)=el ist(2fm) 

A 

0223 

end  i  f 

0224 

c 

i 

0225 

if(nuflag  .eq.  0)  then 

— 

0226 

mhtmx=2 

i; 

0227 

he  1 i st(l)=200 . 

* 

0228 

he  1 i st(2)=0 . 

V 

p-1 
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0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 

0261 

0262 

0263 

0264 

0265 

0266 

0267 

0268 


do  641  n=l,nrspec 

en=a ! og  (coef nu  (n) ) 

cf  I  i st (1 , n)=en+expnu(n) *hc  I  ist(l) 

641  cf  1  i  st(2,  n)=eri4expnu  (n) *hc  I  i  st(2) 
end  i  f 
c 

if(nprof  .eq.  1)  then 

i f  (beta*hpr ime  .eq.  0.  .and.  npath  ,ne.  2)  then 
c  This  is  not  a  presegmented  path,  the  profile  specification 

c  must  be  made  in  the  NAMELIST. 

print  *,’AB0RT  MAIN:  BETA  or  HPRIME  not  input  ’ 
stop 
end  i  f 
nrspec=l 
I htmx=2 

htl i st(l)=200 . 
htl ist(2)=0. 
hprout=hpr i me 
e  I  se 

if(nprof  .eq.  2)  then 

c  Non-exponential  profile,  get  a  value  for  HPRIME 

call  gethpr (wrO, hprout) 
end  i  f 
end  i  f 
c 

c  BEGIN: 
c 

ca I  I  gcpath 
go  to  10 
c 

999  stop 

1000  format(a) 

1001  format(lx, (a)) 

1002  format(8x , i 1) 

1004  format(8x,a) 

1040  f  ormat  ( ’  Input  EIGEN  list:  D=’,f5.2,’  H=’,f5.2/ 
i  *  EIGEN  =’,6(f8.3, ’  ’)/(8x,6(f8 ’  ’))) 

1042  format(’  Reduced  EIGEN  list:  ATNMAX= ’ , f 5 . 1/ 

I  '  EIGEN  =’ , 6 (f 8 . 3,  ’  ’) / (8x,6(f 8 . 3,  ’  ’))) 

end 
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0001  function  cdang(arg) 

0002  complex*16  arg 

0003  real*8  cdang,argr,argi 

0004  argr=drea I (arg) 

0005  argi=dimag(arg) 

0006  cdang=datan2(argi ,argr) 

0007  if(argi  .ge.  O.dO)  return 

0008  cdang=cdang+6 .2831853072d0 

0009  return 

0010  end 


0001 

0002 

0003 

0004 

0018 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 


c 


c 


5 

10 


subroutine  comp  f 
implicit  real+8  (a-h,o-z) 

include  ’common2 . f or ’ 
include  ’common3 . f or ’ 

c=cdcos (theta *zdtr) 
csq=c*c 

s=cds i n (theta*zdtr) 

ssq=s*s 

ca I  I  rbars 

if(rpoly  .eq.  0)  then 
call  i nteg 
else 

ca I  I  uspo I y 
end  i  f 

if(typitr-l)  5,10,15 
f=(rbarll*rll-zone) * (rbar22*r22-zone) 
S  -rbarll*rbar22*rl2*r21 
return 

f=rbarll*rll-7one 

return 

f=rbar22*r22-zone 

return 

end 
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0001 

0002 

0003 

0004 

0020 

0034 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 


subroutine  drvequ 
implicit  real  *8  (a-h,o-z) 
c 

include  ’commonl . for  * 
include  ’common2 . for ’ 
include  ’common3.for ’ 
c 

complex*16  k2i , i I , im, i n.capd, usqd,yud,ysqd, u, usq, 

S  til , t31 , t42, t44 , tl2vrc , tl4 vrc , t32vrc , t34vrc , ct41 , 

$  slla,dlla,sllb,dllb,cl2,c21 , 

S  sl2Jdl2,s21,d21,s22,d22,bll,b22,bl2,b21 

real*8  I sq,msq, nsq, Im, I n ,mn 
real *4  htO 

dimension  cx(3) ,capy (3) ,ysq(3) 

data  dtr/1 . 745329252d-2/ , coef f y /I . 758796dl 1 / , coef f x/3 .  18?357d09/ 
c 

entry  intcmp 

k2 i =dcmp I x (0 . dO , -0 . 5d0*wn) 
s i nd i p=ds i n (cod i p*dtr) 
drcos I =s i nd i p*dcos (az i m*dtr) 
drcosm=s i nd i p*ds i n(az im*dtr) 
drcosn=-dcos  (cod i p*dtr) 
i I =dcmp I x (0 . dO, drcos I ) 
i m=dcmp I x (0 . dO , drcosm) 
i n=dcmp I x (0 . dO, drcosn) 
lsq=drcosl**2 
msq=drcosm**2 
nsq=drcosn**2 
Imsdrcosl *drcosm 
I n=drcosl *drcosn 
mn=drcosm*drcosn 
cO=coeff x/omega**2 
cy=coef fy *magf I  d/omega 
do  1  k=l,nrspec 

cx (k) =cO*charge (k) **2/mrat i o (k) 
capy (k)=cy*charge(k) /mratio(k) 

1  ysq(k)=capy (k) **2 

call  gethpr (100. *wrO,htO) 
topht=htO 
I htmn=l ht 
mhtmn=mht 

if (debug  .le.  1)  return 

print  110 

l=lhtmn 

m=mhtmn 

ht=topht 

10  slopel  =  (ht-htl  i st( Ul))/ (htl  i st ( I ) -ht I  ist(Ul)) 
s I opem=(ht-hc I ist(m>l))/ (he  I ist(m)-hcl ist(m+l)) 
ed=dexp(lnl ist(l+l,l)  +  (lnl i st ( I ,l)-lnl i st( I +1 , 1)) *s I  ope  I ) 
en=dexp(cf I i st(m+l , 1)  +  (cf I i st(m, l)-cf I i st(m+l, 1)) *s lopem) 
capx=ed*cx(l) 
capz=en/omega 
wr=omega*capx/capz 
print  lll^t^d^n^apXjCapZjWr 
i f (ht  .It.  topht)  return 
ht=d 

do  11  j=l , I htmx 
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if(d  .ge.  htlist(j))  go  to  12 

l=j 

do  13  jsmjmhtmx 

if(d  .ge.  hclist(j))  go  to  10 

m=j 

entry  smtrix 
usqd=zero 
yud=zero 
ysqd=zero 

slopel=(ht-htl ist(lht+l))/(htl ist(lht)-htl ist(lht+l)) 
slopem=(ht-hcl ist(mht+l))/ (he  I ist(mht)-hcl ist(mht+l)) 
do  20  k=l,nrspec 
capx=dexp( I n I ist(!ht+l,k) 

♦  (Ini ist(lht,k)-lnl i st ( I ht+1 , k)) *s I  ope  I ) *cx(k) 
capz=dexp  (of  I  i  st(mht-*-l ,  k) 

+(cf I i st(mht,k)-cf I i st (mht+1, k))*slopem) /omega 
u=dcmp I x (1 . dO, -capz) 
usq=u*u 

capd=-capx/ (u* (usq-ysq(k))) 
i f (cdabs(capd)  .gt.  l.d-30)  then 
usqd=usqd+usq*capd 
yud=yud*capy (k) *u*capd 
ysqd=ysqd+ysq(k) *capd 
end  i  f 
conti nue 

crvtrm=a I pha*  (h-ht) 

mll=usqd- I sq*ysqd-crvtrm 

m22=usqd-msq+ysqd-crvtrm 

m33=usqd-nsq*ysqd-crvtrm 

ml2=- i n*yud- lm*ysqd 

m21=  i n*yud- I m*ysqd 

ml3=  im*yud- I n*ysqd 

m31=- im*yud- I n*ysqd 

m23=-i I *yud-mn*ysqd 

m32=  i I *yud-mn*ysqd 

capd=zone/ (zone+m33) 

tll=-s*m31*capd 

tl2vrc=s*m32*capd/c 

tl4vrc=(csq+m33) ♦capd/c 

t31=m23*m31*capd-m21 

t32vrc=c* (m22-m23*m32*capd) /c 

t34vrc=s*m23* capd/c 

ct41=(zone+mll-ml3*m31*capd) *c 

t42=m32*ml3*capd-ml2 

t44=-s*ml3*capd 

slla=tll-»t44 

dlla=tll-t44 

sllb=tl4vrc*ct41 

dllb=tl4vrc-ct41 

sl2=tl2vrc+t42 

dl2=tl2vrc-t42 

s21=t34vrc+t31 

d21=t34vrc-t31 

s22=c+t32vrc 

d22=c-t32vrc 
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0195 

0196 

0197 

0198 

0199 

0200 


i f (ht  ,eq.  topht)  call  i nta I r 
c 

entry  rderiv 
k=0 

do  30  j=l ,7,2 
k=k  +  l 

i f  (dabs( logr (j) )  .gt.  15. dO) 

$  Iogrs(k)=dcmplx(dsign(15.d0, logr(j)) ,0.d0) 

30  rs(k)=cdexp(logrs(k)) 
bll=rll*  (dlla-dllb) 
b22=r22*d22 
bl2=rl2*d21 
b21=r21*sl2 
cl2=rl2*s21 
c21=r21*dl2 
d I lldh=k2i * 

S  (bll+bl2+b21-sllb-sllb+ (rl2*r21*d22+cl2+c21-dlla-dllb) /rll) 

d I22dh=k2i * 

S  (bl2*b21+b22-s22-s22+ (r 12* r21* (dlla-dllb) +bl2+b21+d22) /r22) 

d I 12dh=k2i * 

S  (bll+bl2*b22*slla-sllb-s22+ (rll*sl2*dl2) * (r22+zone) /rl2) 

dl21dh=k2i* 

S  (bll+b21*b22-slla-sllb-s22+ (rll*d21+s21) * (r22+zone) /r21) 

c 

if  (debug  .gt.  2)  then 
print  100, ht,delh, logr, dlrdh 
end  i  f 
return 
c 

100  f ormat (f 9 . 4 , lpel2 . 4 ,4 (lx, 2el2 . 3) /21x,4 (lx, 2el2 .3) ) 

110  format(/’  Electron  density  parameters:  ht  den  nu’, 

S  8x,’x  z  w’) 

111  f ormat(27x, f 7 . 1 , Ip5el0 .2) 
end 


subroutine  extrap 
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0045 
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c 

c  This  routine  sets  up  and  maintains  the  data  sets  for  the  quadratic 
c  extrapolation  of  eigen’s  down  the  propagation  path, 
c 

include  ’commonl .for ’ 
include  ’common2 . for ’ 
c 

logical  brwstr 

complex+8  t(30) ,y (30) , >8(3,30) ,s,tb,stb,capk,coeff , ngsq, 

$  zero/ (0. ,0.)/ , zone/ (1,0.)/, zmp I x i / (0 . , 1 . ) / 

dimension  xs(3) 

equivalence  (el ist,y) , (tl ist,t) 
data  dtr/. 01745329252/ 
c 

if (lx  ,eq.  0)  then 

c  This  is  the  first  point  on  the  propagation  path, 

c  Set  up  constants  and  remove  input  modes  with  attenuation 

c  rates  greater  than  atnmax. 

capk=cmp!x(l . ~.5*alpha*h,0.) 
coeff=cmpl x(0. , 182 .0428*f req) /capk 
c 

c  Get  Brewster  mode 

if (sigma  .It.  l.e-3)  then 
ngsq=cmpl x(epsr,-l . 7975e7*s i gma/f req) 
stb=csqrt(ngsq/ (ngsq+zone) ) *capk 
atten=coef f *stb 
if(atten  . le.  atnmax)  then 

c  The  Brewster  mode  is  contained  within  the  normal  set. 

tb=(90. ,0.) 
brwstr=.fa Ise. 
else 

c  The  Brewster  mode  is  outside  the  normal  set. 

if(atten  .le.  2.*atnmax)  then 
c  The  attenuation  rate  is  not  excessive. 

tb=cmp I x (0 . , -1 . /dtr) *c I og (csqrt (zone-stb**2) ♦zmp I x i *stb) 
brwstr= . true . 
e  I  se 

c  The  attenuation  rate  is  excessive. 

tb=(90. ,0.) 
brwstr= .false, 
end  i  f 
end  i  f 
else 

tb=(90. ,0.) 
brwstr= .false, 
end  i  f 

do  139  k=l , 30 
mode(k)=k 

137  if  (real (y(k))  .gt.  0.)  then 
if (brwstr)  then 

c  If  this  mode  is  near  the  Brewster,  then  keep  it. 

if(abs(  real (y(k)-tb))  .le.  1.  .and. 

$  abs (a i mag(y (k) -tb) )  .le.  .5)  go  to  139 

end  i  f 

atten=coef f *cs i n (y (k) *dtr) 
if(atten  .gt.  atnmax)  then 


EXTRAP 

0086 

do  138  l=k , 30 

0087 

t(!)=t(Ul) 

0088 

138 

y(!)=y(Ul) 

0089 

t(30)=zero 

0090 

y (30)=zero 

0091 

go  to  137 

0092 

end  i  f 

0093 

end  i  f 

0094 

139 

cont i nue 

0095 

return 

0096 

end  i  f 

0097 

c 

0098 

x=rho 

0099 

if(nprint  .gt.  1)  print  1000,x 

0100 

nmds= 1 s 

0101 

do  143  k=l,nmds 

0102 

s=zero 

0103 

do  142  11=1, lx 

0104 

p=l. 

0105 

do  141  12=1, lx 

0106 

if (11  .eq.  12)  go  to  141 

0107 

P=P* (x-xs( 1 2) ) / (xs ( 1 1) -xs( 1 2) ) 

0108 

141 

conti nue 

0109 

142 

s=s+p*ys(ll,k) 

0110 

if(nprint  .gt.  1)  print  1 001, mode (k) ,s 

0111 

t(k)=s 

0112 

143 

y(k)=s 

0113 

c 

0114 

c 

Scan  the  extrapolated  eigen’s  for  invalid  values. 

0115 

do  151  k=l,nmds 

0116 

if (real (y(k))  . le.  0.  .or.  real(y(k))  .ge.  90.  .or. 

0117 

S  aimag(y(k))  .ge.  0.)  then 

print  ERROR  EXTRAP:  Extrapolated  mode’ , mode (k) 

0118 

0119 

lost=l 

0120 

return 

0121 

end  i  f 

0122 

151 

conti nue 

0123 

if(nmds  .It.  30)  then 

0124 

t(nmds+l)=zero 

0125 

y (nmds*l)=zero 

0126 

end  i  f 

0127 

return 

0128 

c 

0129 

entry  xsave 

0130 

c 

0131 

c 

This  entry  point  is  called  after  execution  of  WVGUID. 

0132 

c 

It  updates  the  data  sets  used  to  do  the  quadratic  extrapolation. 

0133 

c 

0134 

x=rho 

0135 

if (lx  .It.  3)  then 

0136 

1 x=l x+1 

0137 

else 

0138 

do  21  1=1,2 

0139 

xs(l)=xs(l*l) 

0140 

do  21  k=l,nmds 

0141 

21 

ys(l  ,k)=ys(Ul,k) 

0142 

end  if 
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0143 

1 s=nmds 

0144 

xs ( 1 x) =x 

0145 

do  25  k=l ,  rwnds 

0146 

25 

ys ( 1 x, k)=y (k) 

0147 

c 

0148 

c 

Keep  all  eigen’s  at  first  input  distance. 

0149 

if (lx  .eq.  1)  then 

0150 

modes=nmds 

0151 

return 

0152 

end  i  f 

0153 

c 

j  is  counter  for  inodes  to  be  output  by  SAVEMC 

0154 

c 

k  is  counter  for  inodes  to  be  used  by  WVGUID 

* 

0155 

j=o 

0156 

k=l 

0157 

251 

if(k  .gt.  nmds)  return 

0158 

j=j*l 

0159 

if(brwstr)  then 

0160 

c 

Check  if  this  mode  is  near  the  Brewster;  if  so,  then  keep  it. 

>  ;\r 

0161 

if(abs(  rea 1 (y (k) -tb) )  .le.  1.  .and. 

0162 

S  abs(aimag(y (k)-tb))  .le.  .5)  go  to  256 

0163 

end  i  f 

0164 

atten=coef f *cs i n (y (k) *dtr) 

0165 

if(atten  .gt.  atnmax)  then 

0166 

c 

Delete  k-th  mode 

w 

0167 

do  253  1=1,4 

A 

0168 

253 

T  term( 1 , j )=zero 

,\V 

0169 

nmds=nmds-l 

0170 

1 s=nmds 

0171 

if (nmds  .eq.  0)  then 

0172 

print  *, ’ERROR  EXTRAP:  All  modes  have  been  deleted’ 

•<  .V 
V.'V 

0173 

lost=l 

0174 

t(l)=zero 

•;  .4 

0175 

y  (l)=zero 

v* 

0176 

return 

1 1 

0177 

end  i  f 

0178 

if(k  .gt.  nmds)  then 

■V 

0179 

257 

modes=modes-l 

i  *  » 

0180 

if (real (T  term(l , modes) )  .ne.  0.)  return 

- 

0181 

go  to  257 

0182 

else 

>> 

0183 

do  254  l=k,nmds 

0184 

mode( 1 )=mode( 1+1) 

'  >*, 

0185 

t(l)=t(Ul) 

■  r 

rt  't* 

0186 

y(l)=y(Ul) 

■S*I 

0187 

do  254  m=l ,  1  x 

0188 

ys(m,  l)=ys(m,  Ul) 

, \ 

0189 

254 

conti nue 

0190 

t(nmds+l)=zero 

0191 

y  (nmds»l)=zero 

T  ' ,.  \ 

0192 

go  to  251 

*■ 

0193 

end  i  f 

,\ul 

i. 

0194 

end  i  f 

*'■ 

0195 

256 

k=k*l 

— 

0196 

go  to  251 

0197 

c 

,  , 

0198 

1000 

format(/’  Extrapolated  EIGEN  list  for  x  =’,f8.3) 

0199 

1001 

format(i5,5x,2f8.3,f 12.3) 

■j 

0200 

end 

,  \ 
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0033 

0034 
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0040 
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subroutine  gcdbr(dl ,c I tl , c I t2, rho, br , i nb) 

Returns  great  circle  distance  and  geographic  bearing  angle 

Input:  DL  is  longitude  of  point  2  minus  longitude  of  point  1 
CLTl  is  co- latitude  of  point  1 
CLT2  is  co- latitude  of  point  2 
INB=0:  RHO  is  computed  and 

BR  from  point  1  thru  point  2  is  computed  at  1 
INB=1 :  RHO  is  input  and 

BR  from  point  1  thru  point  2  is  computed  at  2 

Output:  RHO  is  great  circle  distance  between  the  input  points 

BR  is  geographic  bearing  angle  measured  clockwise  from 
due  North 

All  coordinates,  RHO  and  BR  are  in  radians 
Sign  convention  is  ♦  for  West  and  North 

data  pi /3 . 14159265e0/ , twop i /6 . 28318531eO/ 

reduce(arg)=sign(ami nl (abs(arg) ,1 .) ,arg) 

ccltl=cos(cltl) 
scltl=sin(cltl) 
cclt2=cos(clt2) 
sc  I t2=s i n (c 1 12) 

adl=abs(dl) 
if  (ad I  .ge.  pi)  then 
dl=amod(dl , twopi) 
else 

adl=abs(dl) 
end  i  f 

if(inb  .eq.  1)  then 
if  (rho  .gt.  pi)  then 
gcd=twop i -rho 
else 
gcd=rho 
end  i  f 
end  i  f 

if(abs(cltl)  .le.  l.e-6  .or.  abs(cltl-pi)  .le.  l.e-6)  go  to  10 
if  (ad  I  .le.  l.e-6)  go  to  20 
if (abs(adl-pi)  .le.  l.e-6)  go  to  30 
if  (ad  I  .ge.  pi)  then 
i f (d I  .ge.  0.)  then 
dl=dl-twopi 
else 

d I =d I +twop i 
end  if 
end  i  f 

i f ( i nb  .eq.  0)  then 
cgcd=cc I tl  *cc I t2+sc I tl w sc  I t2*cos (d I ) 
gcd=acos (reduce (cgcd)) 
if (abs(cgcd-l .)  ,le.  l.e-6)  then 
br=0 . 


GC08R 

0058 

br=acos (reduce ( (cc 1 t2-cc 1 tl *cgcd) / (sc  1 tl*s i n (gcd) ) ) ) 

0059 

end  i  f 

0060 

e  1  se 

0061 

if(abs(gcd)  . le.  l.e-6)  then 

0062 

br=0. 

0063 

else 

0064 

br=pi -acos (reduce ((cc 1 tl-cc 1 t2*cos(gcd))/(sc 1 t2*si n(gcd)))) 

0065 

end  i  f 

0066 

end  i  f 

0067 

if(dl  .It.  0.)  br=twopi-br 

0068 

0069 

c 

go  to  40 

0070 

c 

point  1  is  at  one  of  the  poles 

- 

0071 

10 

if(inb  .eq.  0)  gcd=abs (c 1 tl-c 1 t2) 

0072 

if(abs(cltl)  .le.  l.e-6)  then 

0073 

br=pi-dl 

0074 

e  1  se 

0075 

br=d  1 

)  "• 

0076 

end  i  f 

r.1 

0077 

go  to  40 

f  . 

0078 

c 

0079 

c 

coordinates  are  on  same  longitude 

0080 

20 

dc=c 1 tl-c 1 t2 

0081 

i f (dc  . ge .  0. )  then 

■fv 

0082 

br=0 . 

'  4 

0083 

e  1  se 

/l 

0084 

dc=-dc 

» 

0085 

br=pi 

0086 

end  if 

0087 

i f ( i nb  . eq .  0)  gcd=dc 

0088 

go  to  40 

0089 

c 

«* 

0090 

c 

coordinates  are  on  opposite  longitudes 

1. 

0091 

30 

dc=c 1 tl+c 1 t2 

0092 

i f (dc  . 1 e .  pi )  then 

0093 

i f ( i nb  . eq .  0)  then 

0094 

br=0 . 

0095 

else 

« 

0096 

br=p  i 

*' 

0097 

end  i  f 

0098 

e  1  se 

0099 

dc=twop i -dc 

0100 

i f ( i nb  eq.  0)  then 

0101 

br=pi 

0102 

e  1  se 

0103 

br=0 . 

0104 

end  i  f 

0105 

end  i  f 

f  ■ 

0106 

i f ( i nb  .eq.  0)  gcd=dc 

0107 

c 

'Mi 
» 4 

0108 

c 

long  path  calculations 

0109 

40 

i f ( i nb  .eq.  1)  then 

0110 

if(rho  .gt.  pi)  then 

0111 

if (br  .It.  pi)  then 

s 

0112 

br=br+pi 

0113 

0114 

e  1  se 

br=br-p i 
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0115 

end  if 

0116 

end  i  f 

0117 

end  i  f 

0118 

c 

0119 

90 

if(inb  .eq.  0)  rho=gcd 

0120 

return 

0121 

end 
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0001 

0002 

0003 

0004 

0005 

0021 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 


c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


subroutine  gcpath 

sign  convention:  ♦  for  west  and  north,  -  for  east  and  south 

include  ’commonl . for ’ 
include  ’common2.for * 


dimension  profl (50,3) ,prof2(50,3) 
real  I ng, long, I  at, ml ,m2,mgf 
character*72  bcd,preseg 
logical  first 

data  dtr/1 . 745329e-2/ , re/6 . 366/ , a  1 1/80 . / 


mi n=0  -  normal 

=1  -  interpolating  between  preseg  values 

=2  -  last  interpolation  interval 


I  ost=0 
=1 
=2 


no  trouble  with  modes 

dropped  a  mode  in  WVGUID  or  EXTRAP 

all  modes  found  but  one  or  more  changed  significantly 
from  the  extrapolated  values 


nprof=0  -  use  profile  from  MAIN 

1  -  use  exponential  profile 

2  -  read  non-exponential  profiles  along  path 

WARNING:  the  heights  must  match  in  each  profile 


f i rst= . true . 

write(90,*)  ’Failure  at  RHO  1* 


I  x=0 

m  i  n=0 

bta=0 . 

s i g=0 . 

mgf=0. 

sigma 1=0. 

tl ng=tlong*dtr 

tc 1 1= (90 . -t I  at) *dtr 

rhoO=rho 

rhop=rho 

drho=drmi n 

if(rbear  .eq.  720.)  then 

ca I  I  gcdbr((t long- r long) *dtr, tc It, (90.-rlat)*dtr ,gcd, xtr,0) 
brng=xtr/dtr 
i f  ( i gcd  ,eq.  1)  then 
rhomax=gcd*re 
e  I  se 

rhomax=dmax 
end  i  f 
e  I  se 

xtr=rbear*dtr 
brng=rbear 
rhomax=dmax 
end  i  f 


I 


c 

20 

c 


if(npath  eq.  2)  then 
Presegmented  path 
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0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 

0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 

0138 

0139 

0140 

0141 

0142 


read(5,2000,end=900)  preseg 

read (preseg, *)  rho, azm, cdp,mgf , si g,eps , bta , hprm 
if (rho  .eq.  40.)  then 
print  *,’End  of  preseg  data* 
rewind  90 
wr i te (90, 2003) 
go  to  999 

else  if (rho  .gt.  rhomax)  then 

print  *,’DMAX  reached  before  end  of  preseg  data  * 
rew i nd  90 
write (90, 2003) 
go  to  900 
end  i  f 

if  (first)  then 
rho0=rho 
rhop=rho 
end  i  f 

drho=rho-rhop 
i f (drho  .It.  0 . )  then 

print  *, ’ABORT  GCPATH:  Preseg  rhos  out  of  order’ 
go  to  900 
end  i  f 

if(nprof  .eq.  2)  then 
read (3, 2000)  bed 

if (bed (1:8)  .ne.  ’profile  ’  .and. 

S  bed (1:8)  .ne.  ’PROFILE  ’)  then 

print  *, ’ABORT  GCPATH:  PROFILE  control  string  missing’ 
go  to  900 
end  if 

read (bed, 2001)  nn 
if(nrspec  .ne.  max0(l,nn))  then 
print  *, ’ABORT  GCPATH:  Number  of  species  is  incorrect’ 
go  to  900 
end  i  f 

ca I  I  prof i n (3, 1 , 50, npr i nt, nr spec , I htmx , ht I i st, I n I i st) 
i f ( I htmx  .It.  0)  then 
print  *, ’ABORT  GCPATH:  Profile  missing’ 
go  to  900 
end  i  f 

if(lhtmx  .gt.  0)  then 

if (I htmx  .ne.  Ihtmxl)  then 

print  *, ’ABORT  GCPATH:  Number  of  heights  is  incorrect’ 
go  to  900 
end  i  f 

call  gethpr(wrO,hprout) 
end  i  f 
e  I  se 

if(nprof  .eq.  1)  then 
if (bta  .gt.  0.)  then 
beta=bta 
hpr ime=hprm 
end  i  f 

i f (beta*hpr ime  ,eq.  0.)  then 
print  *, ’ABORT  GCPATH:  BETA  or  HPRIME  not  input’ 
go  to  900 
end  i  f 
end  if 
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0143 

c 

Calculate  exponential  profile: 

0144 

Ini i st (1 , l)=cf 1 i st(l , 1) +beta* (ht 1 i st(l) -hpr ime) -9.4517306 

0145 

1 n 1 i st(2, l)=cf 1 i st(2, 1) +beta* (ht 1 i st (2) -hpr ime) -9.4517306 

0146 

hprout=hpr ime 

0147 

end  i  f 

0148 

end  i  f 

0149 

c 

0150 

if(npath  .eq.  1)  then 

‘ 

0151 

c 

Calculate  midpoint  distance: 

0152 

rho= . 5*rhomax 

0153 

else 

0154 

if(rho  .eq.  0.)  then 

0155 

c 

Begin  at  xmtr 

0156 

1 ng=tlng 

0157 

c 1 t=tc  1 1 

0158 

br=xtr 

0159 

end  i  f 

0160 

end  i  f 

* 

0161 

c 

si 

0162 

30 

if(rho  .gt.  0.)  then 

0163 

gcd=rho/re 

0164 

ca 1 1  recvr (t 1 ng, tc 1 t, xtr, gcd ,  1 ng,c 1 1) 

— 

0165 

if(rho  .eq.  20.)  then 

0166 

c 

At  the  antipode  of  the  transmitter. 

»• 

0167 

br=9 . 4247779608-xtr 

0168 

i f (br  .gt.  6.2831853072)  br=br-6. 2831853072 

0169 

else 

’ 

0170 

cal  1  gcdbr(tlng-lng,tclt,clt,gcd,br,l) 

0171 

end  i  f 

0172 

end  if 

t( 

U 

0173 

bpath=br/dtr 

V 

0174 

1 ong= 1 ng/dtr 

h 

t, 

0175 

colat=clt/dtr 

0176 

1 at=90 . -co 1  at 

0177 

c 

0178 

if(sig  .gt.  0.)  then 

0179 

s i gma=s i g 

0180 

epsr=eps 

•K1 

0181 

e  1  se 

'f 

0182 

if(ignd  .eq.  1)  call  ground(long, lat,ncode, sigma, epsr) 

0183 

end  i  f 

0184 

c 

If  conductivity  has  changed,  then  restart  extrapolation 

* 

0185 

if(sigmal  . ne.  sigma  .and.  sigmal  . ne.  0.)  then 

:V 

0186 

mi  n=0 

> 

0187 

1  x=0 

.*: 

0188 

ca 1 1  xsave 

V 

0189 

1  X— 0 

0190 

end  if 

0191 

c 

*'i 

.1 

0192 

i f (mgf  gt.  0.)  then 

K> 

0193 

az i m=azm 

»» 

0194 

cod i p=cdp 

* , 

0195 

magf 1 d=mgf 

r 

0196 

if(magfld  .gt.  l.e-02)  magf 1 d=magf 1 d*l . e-04 

0197 

e  1  se 

•*« 

0198 

cal  1  newmag(0,alt, lng,clt,bmf , d i p, b, br , bp, bt) 

H 

* 

0199 

az im=bpath-bmf /dtr 

;  t,: 
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0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0208 

0209 

0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 


if  (azim  .It.  0.)  then 
az im=az im+360. 

else  if (azim  ,gt.  360.)  then 
az im=az im-360 . 
end  i  f 

cod i p=90 . -d i p/dtr 
magf I d=b*l ,0e-04 
c 

if(mdir  .eq.  1)  then 

c  Reverse  azim  if  contours  for  xmtr  deployment 

az im=az im-180 . 
i f (az im  .It.  0.)  then 
azim=az im+360. 

else  if (azim  .  gt.  360.)  then 
az im=az i m-360 . 
end  i  f 
end  i  f 
c 

end  i  f 

print  1000, rho, I ong, I  at, bpath , az i m, cod i p,magf I d, s i gma , epsr 
c 

40  I ost=0 
x=rho 

ca I  I  extrap 

if (lost  .eq.  1)  go  to  100 
ca I  I  wvgu i d 

iffnmds  .eq.  0  .and.  (rho  .eq.  0.  .or.  npath  .eq.  1  ))  then 
print  *, ’ABORT  GCPATH:  Failure  at  starting  rho’ 
go  to  900 
end  i  f 

if (lost  .eq.  1)  go  to  100 
ca I  I  xsave 

if (lost  .eq.  1)  go  to  100 
rhop=rho 

if  (min  .eq.  1)  go  to  50 
if  (first)  then 

c  Primary  output  file: 

open (un i t= I un i t7 ,status= ’ new’) 

if (year  . eq .  0  .and.  month  .eq.  0  .and.  day  .eq.  0)  then 
wr i te( I un i t7, 1030)  1 1 ong, 1 1  at, brng, beta , hpr ime, path i d 
e  I  se 

wr i te( I un i t7, 1031)  1 1 ong, 1 1  at, brng, beta , hpr ime, 

$  mod (year, 100) , month, day, gmt, path  id 

end  i  f 

f i rst= .false, 
end  i  f 
ca I  I  savemc 
if (npath  . eq .  1)  then 
rew i nd  90 
write (90, 2003) 
go  to  999 
e  I  se 

r* wind  90 

wr . ce (90, 2002)  rho 
end  i  f 

if (rho*. 002  .ge.  rhomax)  go  to  900 
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0257 

rhol=rho 

0258 

al=az im 

0259 

cl=cod i p 

0260 

ml=magf 1 d 

0261 

el=epsr 

0262 

sl=a log (sigma) 

0263 

sigmal=sigma 

0264 

if(nprof  .gt.  0)  then 

0265 

lhtmxl=lhtmx 

0266 

do  48  1=1,1 htmx 

0267 

do  48  m=l,nrspec 

0268 

48 

prof 1(1 ,m)=lnl ist(l  ,m) 

0269 

end  i  f 

0270 

if  (min  .eq.  2)  then 

0271 

mi  n=0 

0272 

go  to  70 

0273 

end  i  f 

0274 

c 

0275 

50 

if  (lost  .eq.  2)  then 

0276 

if (min  .eq.  0)  then 

0277 

drho=amaxl (drho-drmin,drmi n) 

0278 

e  1  se 

0279 

drho=.5* (rho2-rho) 

0280 

end  i  f 

0281 

e  1  se 

0282 

i f (m i n  . eq .  0)  then 

0283 

drho=ami nl (drho+drmin,drmax) 

0284 

else 

0285 

drho=rho2-rho 

0286 

end  i  f 

0287 

end  i  f 

0288 

c 

0289 

70 

if (min  .eq.  0  .and.  npath  .eq.  2) 

go 

to  20 

0290 

rho=rho+drho 

0291 

if(rho+.002  .gt.  rhomax)  then 

0292 

drho=drho-rho+rhomax 

0293 

rho=rhomax 

0294 

end  i  f 

0295 

if (min  .eq.  1)  go  to  120 

0296 

go  to  30 

0297 

c 

0298 

c 

Back  up  on  propagation  path 

0299 

100 

if(rho  .eq.  rhoO)  then 

0300 

print  *, ’ABORT  GCPATH:  Failure 

at 

starting  rho’ 

0301 

go  to  900 

0302 

end  i  f 

0303 

if(min-l)  101,102,103 

0304 

101 

if  (npath  .eq.  2)  go  to  105 

0305 

if(drho  .le.  drmin)  go  to  110 

0306 

nrd=.5*drho/drmi n 

0307 

if(nrd  .eq.  0)  go  to  110 

0308 

drho=nrd*drmi n 

0309 

rho=rhol+drho 

0310 

go  to  30 

0311 

103 

m  i  n=l 

0312 

drho=drhop 

0313 

102 

drho=.5*drho 
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0314 

if(drho  .It.  .015125)  then 

0315 

print  *, 'ABORT  GCPATH:  Backup  interval  is  less  than  0.015125‘ 

0316 

go  to  900 

0317 

end  i  f 

0318 

104 

rho=rhop*drho 

0319 

go  to  120 

0320 

c 

0321 

105 

if  ((rho-rhol)/drmin  .gt.  10.)  then 

0322 

print  *, ’ABORT  GCPATH:  Preseg  interval  too  large  for  efficient 

0323 

{processing' 

0324 

go  to  900 

0325 

end  i  f 

0326 

c 

0327 

c 

Begin  interpolation 

0328 

no 

mi  n=l 

0329 

drho= . 5*  (rho-rhol) 

0330 

if(drho  .It.  .015125)  then 

0331 

print  *, 'ABORT  GCPATH:  Backup  interval  is  less  than  0.015125' 

0332 

go  to  900 

0333 

end  i  f 

0334 

rho2=rho 

0335 

a2=az  im 

0336 

if(a2-al  .gt.  180.)  then 

0337 

a2=a2-360. 

0338 

else  if(a2-al  .It.  -180.)  then 

0339 

a2=a2+360. 

0340 

end  if 

0341 

c2=cod i p 

0342 

m2=magf 1 d 

0343 

e2=epsr 

0344 

s2=a log (sigma) 

0345 

sigma2=sigma 

0346 

if(nprof  .gt.  0)  then 

0347 

1 htmx2= 1 htmx 

0348 

do  111  1=1,1 htmx 

0349 

do  111  m=l,nrspec 

0350 

in 

prof2( 1 ,m)=l n 1 i  st(l ,m) 

0351 

end  i  f 

0352 

rho=rhol*drho 

0353 

C 

0354 

120 

if(rho*.002  .  ge.  rho2)  then 

0355 

c 

End  of  interpolation 

0356 

min=2 

0357 

drhop=drho 

0358 

drho=drmi n 

0359 

rho=rho2 

0360 

az im=a2 

0361 

if (azim  .It.  0.)  then 

0362 

azim=azim+360. 

0363 

else  if (azim  .gt.  360.)  then 

0364 

az im=az im-360 . 

0365 

end  i  f 

0366 

cod i p=c2 

0367 

magf Id=m2 

0368 

epsr=e2 

0369 

s i gma=s i gma2 

0370 

if(nprof  .gt.  0)  then 
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0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 

0386 

0387 

0388 

0389 

0390 

0391 

0392 

0393 

0394 

0395 

0396 

0397 

0398 

0399 

0400 

0401 

0402 

0403 

0404 

0405 

0406 

0407 

0408 

0409 

0410 

0411 

0412 

0413 

0414 

0415 

0416 

0417 

0418 

0419 

0420 

0421 


do  121  1=1,1 htmx 
do  121  m=l,nrspec 

121  Ini i s  t ( I ,m)=prof2(l ,m) 
end  i  f 

el  se 

c  Interpolate 

slope=(rho-rhol) /(rho2-rhol) 
a z im=al  +  s I  ope* (a2-al) 
if (azim  .It.  0.)  then 
az im=az i m+360 . 

else  if (azim  .gt.  360.)  then 
azim=azim-360. 
end  i  f 

cod i p=cl  +  s I  ope* (c2-cl) 
magf I d=ml  +  s I  ope* (m2-ml) 
epsr=el*s I  ope* (e2-el) 
s i gma=exp (sl+s I  ope* (s2-sl)) 
if(nprof  .gt.  0)  then 
do  122  1=1,1 htmx 
do  122  m=l,nrspec 

122  Ini i st ( I ,m)=profl(l ,m) +s I  ope* (prof2( I ,m) -prof 1(1 ,m)) 
end  i  f 

end  i  f 

print  1002,  rho, az i m, cod i p ,magf 1 d , s i gma , epsr 
go  to  40 

c 

900  if(npath  .eq.  2)  then 

903  read (5 , 2000, end=999)  bed 

if (bed (1:5)  .eq.  ’40,0,’)  go  to  999 
go  to  903 
end  i  f 

999  wr i te ( I un i t7 , 1032) 
c I ose (un i t= I un i t7) 

print  ♦,’Execution  terminating  for  this  path’ 
return 

c 

1000  format(/’  Propagation  path  parameters:  rho  long  lat’, 

$  4x, ’bear  azim  codip  magf Id  sigma  epsr’/ 

S  26x, f 10  3,fl0.2,4f9.2,ell . 2, lpell . 2,0pf 8 . 2) 

1002  format(/llx,  ’  Interpolated  path  parameters:  rho  azim’, 

S  5x, ’codip  magf Id  sigma  epsr’/ 

S  38x, f 10  .3, f 10. 2, f9. 2, el 1.2 , lpell .2,0pf 8 .2) 

1030  format(’sw  xmtr ’ , f 7 . 1 ,2f 6 . 1 , ’  prof ’ ,f5.2,f5.1/a80) 

1031  format(’sw  xmtr ’ , f 7 . 1 , 2f 6 . 1 , ’  prof ’ , f5 .2,f5 . 1 , 

$  ’  (\3(i2.2,7’),f4.1,’)7a80) 

1032  format(’r  40.’) 

2000  format(a72) 

2001  f orriat  (8x ,  i  1) 

2002  format(f6.3) 

2003  f ormat ( '40  ’) 
end 


subroutine  gethpr (wr,hpr) 


Routine  to  determine  the  height  where  omega  sub  r  is  a 
specific  value.  The  value  returned  is  to  nearest  km. 


include  ’commonl .for’ 


data  coeffx/3.182357e9/ 


Start  at  the  bottom  of  the  profile  and  work  up. 


I ht= ! htmx-1 
mht=mhtmx-l 

ht=am i nl (ht I i st ( I htmx) , he  I i st (mhtmx)  ) 
i f ( I ht  .gt.  1  .and.  ht  .ge.  htl ist(lht-l))  then 


I ht= I ht— 1 
go  to  10 
end  if 

if(mht  .gt.  1  .and.  ht  .ge.  he  I ist(mht-l))  then 


mht=mht-l 
go  to  12 
end  i  f 

slope  I = (ht-ht I ist(lht+l))/ (htl i st ( I ht) -ht I ist(lht+l)) 
slope  m=(ht-hcl  ist(mht+l))/ (he I  ist(mht)-hcl  ist(mht«-l)) 
sum=0. 

do  14  n=l,nrspec 

dn=exp(lnl i st( I ht+1 ,n)  +  ( I n I i st( I ht, n) -  I n I i st( I ht+1 , n)) *s lope  I) 
cf=exp(cf I ist(mht+l,n)+(cf I i st(mht, n)-cf I ist(mht+l,n))*slope  m) 
sum=sum+coeffx*dn/(mratio(n)*cf) 
if (sum  .gt.  wr)  then 
hpr=ht 

if(nprint  .gt.  1)  then 
print  * 

print  *,’GETHPR:  wr,  hpr=',wr,hpr 
end  if 
return 
end  i  f 
ht=ht+l . 
go  to  10 
end 


subrout i ne  ground (x long, x I  at, ncode, s i gma , epsr) 

Returns  conductivity  code,  conductivity,  dielectric  constant 


0001 

0002  c 

0003  c 

0004  c 

0005  c 

0006  c 

0007  c 

0008  c 

0009  c 

0010  c 

0011  c 

0012  c 

0013  c 

0014  c 

0015  c 

0016  c 

0017  c 

0018 
0019  1 

0020  1 
0021  1 
0022  c 

0023 
0024 
0025 
0026 
0027 
0028 
0029 
0030 

0031  1 

0032 
0033 
0034 

0035  2 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046  11 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 


Input:  XLONG  is  West  longitude  in  degrees 

(i.e.,  -117.3  for  117  degrees,  18  minutes  East) 

XLAT  is  North  latitude  in  degrees 

(i.e.,  32.8  for  32  degrees,  48  minutes  North) 

Output:  NCODE  is  conductivity  code  from  GRNDMAP.DAT 

(ncode=0  is  sea  water,  =1  is  ice;  see  DATA  below) 

SIGMA  is  mho/m 

EPSR  is  the  dielectric  constant 

A  I i st  of  sigma  and  epsr  is  also  placed  into  a  common. 
Requires:  GRNDMAP.DAT 

i nc I ude  * [305021 . jaf I i b] data_f i I es .f or/ list’ 
character*40  grndSd/ 'userSd i skS3 : [305021 .jaf I i bjgrndmap . dat ’/ 
character *40  i tsnSd/ ’ userSd i sk$3 : [305021 . jaf I i b] i tsnoi se .dat’/ 
character *40  wr I dSd/ ’userSd i sk$3 : [305021 .jaf I i b] wor I d .dat ’/ 

common/grnd$/sss(10) , rrr (10) 
logical  f i rst/ . true . / 

dimension  !code(361) ,map(4530) ,ss(10) ,rr(10) 

data  ss/1 . e-5,3.e-5,l .e-4,3.e-4,l .e-3,3.e-3,l .e-2,3.e-2, .1,4./ 

data  rr/5., 5. ,10. ,10. ,15. ,15. ,15. ,15. ,80. ,81./ 

if (first)  then 

open (un i t=8 , f i I e=gr  ndSd , sta tus= *  o I d ’ , readon I y) 
read  (8, 1)  lcode,map 
format(9i8) 
close(unit=8) 
do  2  1=1,10 
sss(l)=ss(l) 
rrr(l)=rr(l) 
f i rst= .false, 
end  i  f 
ph i =x I ong 

if  (phi  .  gt.  180.)  then 
ph  i  -ph  i  -360 . 
e  I  se 

if  (phi  .It.  -180.)  ph  i  =ph  i  -*-360 . 
end  i  f 

if(abs(phi)  gt.  180.  .or.  abs(xlat)  .gt.  90.01)  then 
print  11 , x I ong, x I  at 

format(/’  **********  Error  in  GROUND:  Xiong  Xlat' 

S  /26x, 2f 9 . 2) 

stop 
end  i  f 

lat=181 . -2. *xlat 
if (I  at  gt.  360)  lat=360 
I ong=361 . -2 . *ph i 
if (long  .gt.  720)  long=l 
1 1= I code( I  at) 

1 2= I  code ( I at+1) -1 
do  21  1=11,12 
maplml=map(l)/10000 
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GROUND 


0058 

map  1 =map( 1 ) -10000*map tml 

0059 

mlong=mapl /10 

0060 

if(mlong  . ge.  long)  go  to  31 

0061 

21 

continue 

0062 

31 

ncd=map 1 -m 1 ong* 10 

0063 

mlong=maplml/10 

0064 

if(mlong  .It.  long)  go  to  41 

0065 

ncd=map 1 ml-m 1 ong*10 

0066 

41 

if(ncd  .It.  0  .or.  ncd  .gt.  9)  then 

0067 

print  51,xlong,xlat,  long,  lat,  11, 12,  ml  ong,  ncode 

0068 

51 

format(/'  **********  Error  in  GROUND:  Xiong  Xlat', 

S  *  long  lat  11  12  ml ong  ncode' 

0069 

0070 

S  /26x , 2f 9 . 2 , 4 i 6 , 2 i 8) 

0071 

stop 

0072 

end  i  f 

0073 

ncode=ncd 

0074 

if (ncode  .eq.  0)  ncode=10 

0075 

sigma=ss (ncode) 

0076 

epsr=  rr (ncode) 

0077 

return 

0078 

end 

c 


subroutine  intalr 
implicit  real  *8  (a-h,o-z) 


0001 

0002 

0003 

0004 

0018 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 

0056 

0057 

0058 
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0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 


c 

10 


c 


30 


40 


50 


60 


70 


90 


include  ’common2 . for ’ 
include  *common3 . f or ’ 

comp  I  ex* 16  q,p, t,dll , dl3, d31 , d33, del ta , f  nsq, f root, 

S  coml ,com3,csqm22,csqm33,b3,b2,bl,bO 

dimension  phase!  (8),  phase2(8),  p(2),  t(2),  q(4) 
equivalence  ( I ogrll , phasel  (1) ) 
data  pi/3. 141 592653d0/ , teop i /6 . 283185307d0/ 

if(isotrp-l)  10,100,102 

coml=zone*mll 

com3=zone+m33 

csqm22=csq+m22 

csqm33=csq+m33 

b3=0 . 25d0*s* (ml3*m31) /com3 

b2= (-csqm33*coml*ml3*m3l-com3*csqm22*m23*m32) / (6 . d0*com3) 
bl=s* (ml2*m23*m32*m21-csqm22* (ml3*m31) ) / (4 . d0*com3) 
b0= (coml*csqm22*csqm33+ml2*m23*m31*m32*m21*ml3-ml3*m31*csqm22 
S  -coml*m23*m32-ml2*m21*csqm33) /com3 

cal  I  qartic(q,b3,b2,bl,b0, debug, newq) 

do  30  n=l , 2 

dll=zone*mll-q(n) **2 

dl3=ml3*s*q(n) 

d31=m31+s*q(n) 

d33=zone+m33-s**2 

del ta=dll*d33-dl3*d31 

p (n) = (-ml2*d33+dl3*m32) /de I ta 

t (n) =q(n) *p(n) -s* (~dll*m32*ml2*d31) /del ta 

pyntng=t(n)*dconjg(p(n))*q(n) 

if(pyntng  .It.  0.)  print  201 , theta ,q(n) ,pyntng 
conti nue 

de I ta= (t(l)*c+p(l))*(c*q(2))-(t(2) *c*p(2)) * (c*q(l)) 
rll  =((t(l)*c-p(l))*(c*q(2))-(t(2)*c-p(2))*(c*q(l)))/delta 
r22  =((t(l)*c+p(l))*(c-q(2))-(t(2) *c+p(2) ) * (c-q(l)) ) /del ta 
rl2  =-2 .dO*c* (t(l) *p(2)-t(2) *p(l))/del ta 
r21  =-2.d0*c*(q(l)-q (2) ) /de I ta 
logrll=cdlog(rll) 
logrl2=cdlog(rl2) 
logr21=cdlog(r21) 
logr22=cdlog(r22) 
i f (ad j f I g  .eq.  1)  then 
do  70  n=2,8,2 

i f (phasel (n) -phase2(n)  .le.  pi)  go  to  60 
phasel (n)=phasel (n) -twop i 
go  to  50 

i f (phase2 (n) -phasel (n)  .le.  pi)  go  to  70 
phasel (n) =phasel (n) *twop i 
go  to  60 
conti nue 
end  i  f 

do  90  n=2, 8, 2 
phase2(n)=phasel (n) 
if  (debug  .gt.  2)  print  202 
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return 

c 

100  i r=l 

f nsq=zone*mll 
f root=cdsqrt(csq+mll) 
go  to  106 

101  rll= (f nsq*c-f root) / (f nsq*c  +  f  root) 
r 22=(c-f root) / (c*f root) 

go  to  105 
c 

102  i r=2 

f nsq=zone*mll 

f root=cdsqrt(csq+mll+ml3**2/f nsq) 
go  to  106 

103  coml=(s*f  root+ml3)/(s*fnsq+ml3) 
rl 1= (c-coml) / (c+coml) 

i  r=3 

f root=cdsqr t (csq*m22) 
go  to  106 

104  r22=(c-f root) /(c+f root) 

105  rl2=(l .d-20,0 . dO) 
r21=(l .d-20,0.d0) 
go  to  40 

c 

106  if (di mag (f root)  .gt.  O.dO)  froot=-froot 
i f  (i r-2)  101,103,104 

c 

201  format(’  for  thetas’ , f 7 .4, f9 .4, *  q=’ ,lp2ell .3, 

S  ’  poynti ng(z)=’ ,ell .3) 

202  format(/4x, ’ht’ , 7x, ’del h ’) 
end 
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0001 

subroutine  integ 

0002 

implicit  real  *8  (a-h,o-z) 

0003 

c 

0004 

include  ’commonl . f or  * 

0020 

include  ’common2 .for ’ 

0034 

include  ’comroon3 . f or * 

0059 

c 

0060 

real *8  logrO 

0061 

integer  sflag 

0062 

dimension  logrO (8),  dtrdh0(8),  dlogr0(8),  dlogrl(8),  dlogr2(8) 

0063 

data  dlhmin/1  . 953l25d-3/,dlgrmx/l.d20/ 

0064 

c 

0065 

factor=10.d0** (-prec) 

0066 

emax=f actor *3 .  dO 

0067 

emi n=f actor* . 3d0 

0068 

ht=topht 

0069 

1 ht=lhtmn 

0070 

mht=mhtmn 

0071 

delh=3. 125d-2 

0072 

svdel h=del h 

0073 

if(debug  .gt.  2)  print  200, theta 

0074 

ca 1 1  smtr i x 

0075 

c 

0076 

c 

runge  kutta 

0077 

10 

sf 1 ag=0 

0078 

if (debug  .gt.  2)  print  201 

0079 

11 

if(lht  .It.  Ithmx  .and.  ht  .le.  htl i st( 1 ht*l) )  then 

0080 

lht=lht«-l 

0081 

go  to  11 

0082 

end  i  f 

0083 

13 

i f (mht  .It.  mthmx  .and.  ht  .le.  hcl ist(mht«-l))  then 

0084 

mht=mht+l 

0085 

go  to  13 

0086 

end  i  f 

0087 

if(ht-delh  .It.  htl ist(!ht+l))  then 

0088 

sf lag=l 

0089 

saveht=htl ist(lht*l) 

0090 

del h=ht-saveht 

0091 

end  i  f 

0092 

if(ht-delh  .It.  d)  then 

0093 

sf 1 ag=l 

0094 

saveht=d 

0095 

del h=ht-saveht 

0096 

end  i  f 

0097 

do  30  i=l,8 

0098 

Iogr0(i)  =  !ogr  (i) 

0099 

30 

dl rdhO(i)=dl rdh(i) 

0100 

c 

0101 

c 

Try  again 

0102 

40 

do  50  i=l#8 

0103 

d 1 ogrO(i )=-d 1 rdhO(i ) *de(  h 

0104 

50 

1 °9r  ( i )  =  1 ogrO ( i ) *0 . 5d0*d 1 ogrO ( i ) 

0105 

ht=ht-0.5d0*delh 

0106 

ca 1 1  smtr i x 

0107 

do  60  i =1 , 8 

0108 

d 1 rdh (i )=ds i gn (dmi nl (d 1 grmx,dabs  (d 1 rdh(i))) ,dl rdh(i)) 

0109 

d logrl (i)=-d 1 rdh(i) *delh 

a 

10 

60 

1 ogr ( i )= 1 ogr0( i ) *0 . 5d0*d 1 ogrl ( i ) 

01 

u 

ca 1  I  rder i v 

01 

12 

do  70  i=l ,8 

2 

13 

d 1 rdh ( i )=ds i gn (dmi nl (dlgrmx,dabs (d 1 rdh(i ))) ,d 1 rdh( i )) 

01 

14 

dlogr2(i)=-dl rdh(i)*delh 

01 

15 

70 

logr(i)  =  logr0(i)-»dlogr2(i) 

01 

16 

ht=ht-0.5d0*del h 

17 

ca 1 1  smtr i x 

1 1 

18 

error=0.d0 

19 

do  80  i=l,8 

d 1 rdh ( i ) =ds i gn (dm i nl (d 1 grmx , dabs (d 1 r dh ( i ) ) ) , d 1 rdh ( i )  ) 

01 

dlogr4=((-dl rdh ( i ) *del h*d 1 ogr0( i )) /2 .dO+d 1 ogrl (i) +d 1 ogr2(i ))/3.d0 

2 

1 ogr ( i ) = 1 ogrO ( i ) +d 1 ogr4 

01 

80 

error=error+ (d 1 ogr2(i)-dlogr4) **2 

01 

24 

error=dsqrt(error/8.d0) 

2 

if (error  .It.  emax  .or.  delh  . le.  dlhmin)  go  to  100 

01 

26 

sf 1 ag=0 

01 

ht=ht+del h 

28 

delh=0.5d0*delh 

!  1 

29 

if (delh  .It.  dlhmin)  delh=dlhmin 

go  to  40 

100 

ca 1 1  rder i v 

01 

if (error  .It.  emin)  delh=2.*delh 

2 

if(sflag  .eq.  1)  then 

01 

34 

delh=svdelh 

2 

htssaveht 

01 

36 

end  i  f 

EH 

svdel h=del h 

§ 

38 

if(ht  .gt.  d)  go  to  10 

01 

39 

return 

eh 

40 

c 

01 

41 

200 

format (/’  DEBUG:  theta  =*,2f9.4) 
format(’  ’) 

01 

42 

201 

3! 

43 

end 
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subroutine  iterat 
implicit  real  *8  (a-h,o-z) 


This  routine  drives  the  iteration  to  obtain  solutions  to  the 
modal  equation. 


include  ’common2.for ’ 
include  ’common3 . f or ’ 


complex*16  thetaO, fO, d I thta 
real*4  absr,absi 


nr i ter=0 

if (debug  .gt.  1)  then 
if(rpoly  .eq.  0)  then 
print  300 


print  301 
end  i  f 
print  302 
end  i  f 

Store  the  starting  angle 

thetaO=theta 

theta=theta-dthta 

ca I  I  comp  f 

fO=f 

theta=theta+dthta 
ca I  I  comp  f 

Store  the  magnitude  of  the  f-f unction  for  the  starting  angle 
if (nr iter  .eq.  0)  fmagO=cdabs(f ) 
nr i ter=nr i ter+1 
dfdtht=(f-fO) /dthta 
d I thta=-f /df dtht 
if (debug  .gt.  1)  then 
fmag=cdabs(f) 

print  303, theta , fmag, d I thta , df dtht 
end  i  f 


absr=dabs(drea I (dlthta)) 
abs i =dabs (d imag(d I thta) ) 

if(absr  .gt.  thtinc)  d I thta=d I thta* (tht i nc/absr) 
if(absi  .gt.  thtinc)  d I thta=d I thta* (tht i nc/abs i ) 
theta=theta+d I thta 
if (nr  iter  .It.  maxitr  .and. 

S  (absr  .gt.  lub(l)  .or.  absi  .gt.  Iub(2)))  go  to  10 


nr i ter=nr i ter+1 
fO=f 

ca  I  I  comp  f 
df dtht=(f-f0) /d I thta 
if (debug  .gt.  1)  then 
fmag=cdabs(f) 
d I thta=zero 

print  303, theta , fmag, d I thta , df  dtht 
end  i  f 

if(rpoly  .eq.  1)  then 

Test  the  magnitude  of  the  f-function  of  the  final  angle 
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fmag=cdabs(f) 
if(fmag  .gt.  fmagO)  then 
print  304, fmagO ,f mag 
theta=thetaO 
end  i  f 
end  i  f 

if(typitr  ,gt.  0)  then 
if(typitr  .eq.  1)  then 
dfdtht=(rbar22*r22-zone) *dfdtht 
else 

df dtht=(rbarll*rll-zone) *df dtht 
end  i  f 
end  i  f 
return 

f ormat( ’Olterati ons :  exact’) 

format( ’Olterati ons :  inexact') 

f ormat(8x, ’ rea I  imag  f  mag  d  real  d  imag’,5x, 

$  ’dfdt  real  dfdt  imag’) 

format(5x,2f8 .4 ,lpel2 .3,2(lx,2ell .3)) 

f ormat( ’  Warning  HERAT:  During  RP0LY=1,  starting  fmag  (’, 

S  lpelO.4,’)  is  smaller  than  final  fmag  (’ ,lpel0.4, ’) ’) 
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subroutine  mdhnkl  (z , hi , h2, hlprme, h2prme, theta , i dbg) 


implicit  complex*16  (a-h,o-z) 
complex*16  i ,mpower ,mterm 
rea 1*8  a, b,c,d, cap , parti , part2, zmag 
character*4  idbg 

dimension  a (30) ,  b(30),  c(30),  d(30),  cap (30),  parti (2),  part2(2) 
equivalence  (parti , term4) ,  (part2,sum4) 


data  a 

$ 

S 

$ 

s 

s 

$ 

s 

$ 

$ 

$ 

s 

s 

s 

s 

data  b 

$ 

S 

$ 

s 

s 

s 

s 

s 

$ 

s 

s 

s 

$ 

s 

data  c 

S 

$ 

s 

s 

$ 

s 

s 

s 

$ 

$ 

s 

s 

s 

s 

data  d 

S 

s 

s 


/  9 . 3043671692922944819d-01 

2 .06763714873!S209897d+02 
8 . 7021 7655 1 900761 7234d*02 
5 . 4168543740434246542d*02 
9 . 345849506631 1674231d+01 
6 . 1210004300561072794d+00 
1 ,8401275944132116616d-01 
2 . 8842080097260218300d-03 
2 . 5827494893312753646d-05 
1 ,4155736366074870734d-07 
5 . 0110220346327933889d-10 
1 . 1961806496091228666d-12 
1 ,9948392989517716388d-15 
2 . 3938095525516785112d-18 
2 . 1194208514407528762d-21 
/  6 . 7829872514427588456d-01 

5 . 38332321 54 307 609704d+01 
1 .5337103177865415841d*02 
7 . 4742218215718400631d*01 
1 . 0785312873841039006d*01 
6 . 1360373635097223595d-01 
1 .6422939954686564465d-02 
2 . 3316778764072130571d-04 
1 . 9156708045016374595d-06 
9 . 72861 24416697769730d-09 
3 . 2160808603234314644d-ll 
7 . 2151886229105003778d-14 
1 . 1 368553061 173507104d-16 
1  2946984700995355913d-19 
1 . 0920223904914870636d-22 
/  4 . 6521835846461472410d-01 

2 . 584 54643591 45262382d*01 
6 .  2158403942148298012d*01 
2.70842718702171232280*01 
3 . 5945575025504490022d*00 
1  9128126343925335199d-01 
4 . 8424410379295043444d-03 
6 . 55501 82039227768583d-05 
5 . 16549897866255071 19d-07 
2 . 52 781 00653705 126277d- 09 
8 . 082293604 24644091 57d- 12 
1 ,7590891906016512675d-14 
2 . 695728782367258964 Id- 17 
2 . 99226 19406895981 31 5d-20 
2 . 464 4 4285051 25033375d-23 
/  6 . 7829872514427588456d-01 

3 . 768326250801 5326776d+02 
1 . 9938234131225040548d+03 
1 .4201021460986496090d+03 


3 . 1014557230974314911d+01 , 
5 . 7434365242545027449d*02, 
8 . 2877871922864397320d*02, 
2 . 57945446383020221 lld*02, 
2 . 6626351870744066662d*01 , 
1 . 1592803844803233472d*00, 
2.4833030963741048003d-02, 
2 . 9133414239656786138d-04, 
2 . 0256858739853140063d-06 , 
8 . 8695090013000443124d-09, 
2 . 5658074934115685526d-ll , 
5 . 0988092481207283185d-14 , 
7 . 1886100863l26905797d-17 , 
7 ,3883010881224645255d-20, 
5 . 6653858632471341093d-23/ 
1 . 1304978752404598033d+01 , 
1 . 1962940478735024376d+02 , 
1 ,2780919314887846509d*02, 
3 . 23559386215231 17060d+01 , 
2 . 8532573740320209005d+00 , 
1 ,0937678009821251966d-01, 
2 . 1 05505122395713391 ld-03, 
2 . 2528288660939256561d-05 , 
1 ,4446989475879618839d-07, 
5 . 8854279743918795891d-10, 
1 ,5952782045255116351d-12, 
2 . 9876557444763976717d-15, 
3 . 9889659863766691603d-18, 
3 . 8985199340546088228d-21 , 
2 . 8527230681595795812d-24/ 
6 ,2029114461948629822d+00, 
5 . 2213059311404570392d*01 , 
4 . 8751689366390821897d+01 , 
1 . 1215019407957400909d*01 , 
9 . 1815006450841609147d-01 , 
3 . 3122296699437809740d-02, 
6 . 0568368204246458321d-04 , 
6 . 1985987743950608612d-06, 
3. 82204881 884021 50986d-08, 
1 . 5033066103898380141d-10, 
3 . 9473961 437101054471d-13, 
7 . 1 81421 4762263778920d-16, 
9 . 3358572549515461865d-19, 
8 . 901 567576051 1620701d-22 , 
6 . 3656020935361057409d-25/ 
4 . 521991 5009618392131d+01 , 
1 . 1962940478735024344d+03, 
2  0449470903820554375d*03, 
7 . 1 1 83064967350857463d*02 , 
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0058 

0059 

0060 

0061 

0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 


c 


50 

60 


S 

s 

$ 

s 

$ 

s 

s 

s 

$ 

s 

s 

data 

S 

$ 

$ 

s 

$ 

$ 

s 

$ 

s 

s 

s 

s 

s 

s 

data 

data 

data 

data 

data 

data 

data 

data 


2 . 6963282184602597492d+02, 
1 . 9021 71 58268801 39294d +01 , 
6 . 0764877832340288572d-01 , 
1 ,0026214868551016149d-02, 
9 . 3867869420580235442d-05, 
5 . 3507368429183773360d-07 , 
1 ,9618093247972931935d-09, 
4 . 8341763773500352579d-12, 
8 . 2990437346566602039d-15, 
1 ,0228117913786331176d-17, 
9 . 2821903191776400453d-21 , 
cap  /  1.0416666666666666663d-01, 

1 ,2822657455632716019d-01; 
8 . 8162726744375764874d-01 , 
1 ,4995762986862554546d+01 , 
4 . 74451 53886826431 887d+02, 
2 . 4086549640874004 605d+04 , 
1 . 7919020077753438063d+06, 
1.83707379676330729788+08, 
2 . 482751 9375935888472d+10, 
4.2771126865134715582d+12, 
9 . 1486942234356396792d+14 , 
2 . 3788844395175757942d+17 , 
7 . 3900049415704853993d+19, 
2 . 7030825930275761623d+22 , 
1 . 1498937014386333524d+25, 
i / (0 . dO , 1 . dO) / 


7.9891 20647289o.^o* iiu>G, , 
3 . 7188105233392256682d+00, 
8 . 4220204895828535644d-02 , 
1  0363012784032058021d-03, 
7 . 5124345274574017960d-06, 
3 . 4135482251472901638d-08, 
1 .0209780508963274472d-10, 
2 . 09135902113347837238-13, 
3 ,0316141496462685641d-16, 
3 . 1967863459247792364d-19, 
2 . 5103962999804300309d -22 / 
8 . 3550347222222222116d-02 , 
2 . 9 1 8490264641 404631 5d -01 ; 
3 . 3214082818627675264d+00, 
7 . 8923013011586517530d+01 , 
3.2074  9009089066 1 9004d +03 , 
1 ,9892311916950979121d+05, 
1 . 7484377l80034121023d+07, 
2 . 0679040329451 551 508d+09 , 
3 . 1669454981734887315d+ll , 
6 . 0971132411392560749d+13, 
1 . 4413525170009350101d+16, 
4 . 1046081600946921885d+18, 
1 . 38592200046039431 41d+21 , 
5 ,4747478619645573335d+23, 
2 . 5014180692753603969d+26/ 


one/ (1 ,d0,0 .dO)/, two/ (2. d0,0 . dO) / , zero/ (0 . d0,0. dO) / 
root3/ (1 . 732050807 56888d0 , 0 . dO) / 
a  I pha/ (8 . 53667218838951d-l ,0 . dO) / 

constl / (  2 . 588l9045102522d-01 , -9 . 65925826289067d-01) / 
const2/(  2 . 58819045102522d-01 ,  9 ,65925826289067d-01) / 
const3/ (-9 . 65925826289067d-01 ,  2 . 58819045102522d-01) / 
const4/ (-9 . 65925826289067d-01 , -2 . 58819045102522d-01) / 


zpower=one 
sum3=zero 
sum4=zero 
zmag=cdabs  (z) 

if(zmag  .gt.  6.1d0)  go  to  70 

suml=zero 

sum2=zero 

zterm=-z+*3/ (200 . d0,0 . dO) 
do  50  m=l , 30 

suml=suml +dcmp I x (a (m) , 0 . dO) *zpower 
sum2=sum2+dcmpl x (b(m) , 0 . dO) *zpower 
sum3=sum3+dcmp I x (c (m) , 0 . dO) »zpower 
ternn4=dcmp  I  x  (d  (m)  ,  0 .  dO)  *zpower 
sum4=sum4+term4 

i f (dabs (parti  (1) )  .le.  1  d-17*dabs(part2(l))  .and. 

$  dabs (parti (2) )  .le.  1 .d-17*dabs(part2(2)))  go  to  60 
zpower=zpower*zterm 
gm2f=  i  *  (z»sum2-two*SLiml)  /root3 
gpmfp=i + (sum4+two*z*z*sum3)/root3 
hl=z*sum2+gm2f 
h2=hl-two*gm2f 
hlprme=sum4+gpmf  p 
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0115  h2prme=hlprme-two*gpmf p 

0116  go  to  999 

0117  70  mpower=one 

0118  suml=one 

0119  sum2=one 

0120  rtz=cdsqrt(z) 

0121  sqrtzb=rtz*z 

0122  zterm= i /sqrtzb 

0123  mterm=-zterm 

0124  dm=zero 

0125  term3=one 

0126  do  80  m=l , 30 

0127  zpower=zpower*zterm 

0128  mpower=mpower*mterm 

0129  dm=dm+one 

0130  terml=dcmp I x(cap(m) ,0 . dO) *zpo#er 

0131  term2=dcmp I x (cap (m) ,0 . dO) *mpower 

0132  i f  (cdabs  (term2/term3)  .ge.  l.dO)  go  to  81 

0133  suml=suml*terml 

0134  sum2=sum2+term2 

0135  sum3=sum3+dm*terml 

0136  term4=dm*term2 

0137  sum4=sum4+term4 

0138  i f  (dabs  (parti (1) /part2(l) )  . le.  l.d-17  .and. 

0139  $  dabs  (parti (2) /part2 (2) )  . le.  l.d-17)  go  to  81 

0140  80  term3=term2 

0141  81  zterm=(-1.5d0,0.d0)/z 

0142  sum3=sum3*zterm 

0143  sum4=sum4*zterm 

0144  terml=  ( (-0 . 25d0 , 0 .  dO)  -  i  *sqr  tzb)  /z 

0145  term2= ( (-0 . 25d0 , 0 . dO) ♦ i *sqr tzb) / z 

0146  expl=cdexp ( (0 . dO , 0 . 666666666666666667d0) *sqrtzb) 

0147  exp2=constl*expl 

0148  exp3=const2/expl 

0149  exp4=const3*expl 

0150  exp5=const4/expl 

0151  zterm=a I pha/cdsqrt(rtz) 

0152  term4=z 

0153  if  (parti (1)  ,ge.  O.dO  .or.  partl(2)  .ge.  O.dO)  go  to  90 

0154  hl=zterm* (exp2*sum2+exp5*suml) 

0155  hlprme=zterm* (exp2* (sum2*term2+sum4) +exp5* (suml*terml*sum3) ) 

0156  go  to  110 

0157  90  hl=zterm*exp2*sum2 

0158  hlprme=zterm*exp2* (sum2*term2+sum4) 

0159  110  if(partl(l)  .ge.  O.dO  .or.  partl(2)  .It.  O.dO)  go  to  120 

0160  h2=zterm* (exp3*suml+exp4*sum2) 

0161  h2prme=zterm* (exp3* (suml*terml+sum3) +exp4* (sum2*term2*sum4) ) 

0162  go  to  999 

0163  120  h2=zterm*exp3*suml 

0164  h2prme=zterm*exp3* (suml*terml+sum3) 

0165  c  calculate  wronskian  as  partial  check  on  validity 

0166  999  sum4=hl*h2prme-hlprme*h2 

0167  i f  (dabs  (part2 (1) )  .le.  l.d-8  and. 

0168  S  dabs(part2(2)+l . 4 5 74 954 41 04046 ldO)  .le.  l.d-8)  go  to  1000 

0169  print  1001 , sum4 , theta , i dbg 

0170  1000  return 

0171  1001  format(’  *****  possible  error  in  mdhnk I :  w  =  ’,lp2el5.6, 
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0172  S  ’  for  theta  =  \0p2fl0.4,’  at  \a4) 

0173  end 
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0001  subroutine  newmag ( j , r , phi j , thet, bmf ,di p, b,br , bp, bt) 

0002  c 

0003  c  Returns  parameters  of  the  geomagnetic  field 

0004  c 

0005  c  Input:  J=0:  Use  spherical  earth 

0006  c  J=1 :  Use  spheroidal  earth 

0007  c  R  is  altitude  in  km 

0008  c  PHIJ  is  West  longitude  in  radians 

0009  c  THET  is  co- latitude  in  radians 

0010  c 

0011  c  Output 

0012  c 

0013  c 

0014  c 

0015  c 

0016  c 

0017  c 

0018  dimension  g(10,10),  bm(10) 

0019  da ta  g/ . 0, 3 . 032193e04 , 2 . 522093e03 , -3 . 285459e03 , -4 . 170639e03 , 1 . 6928 

0020  $19e03 , -6 . 684202e02 , -1 . 900312e03 , -2 . 405232e02 , -9 . 358495e02 , -5 . 75507 

002 1  $0e03 , 2 . 1 31 549e03 , - 5 . 206994e03 , 6 . 237642e03 , -4 . 496227 e03 , -3 . 650850e0 

0022  S3 , - 1 . 241578e03 , 2 . 029996e03 , -4 . 463745e02 , -3 . 659410e02, 3 . 495705e03, - 

0023  S 1 . 085898e02 , - 1 . 369823e03 , -2 . 514676e03 , -1 . 943789e03 , -1 . 836598e03 , -1 

0024  S . 31 3045e02 ,  -1 . 626874e02 , 4 . 917246e02 , -8 . 068787e02 , 1 . 220352e03 ,-4.75 

0025  S3 1 92e02 , 1 . 392784e02 , -6 . 836385e02 , 8 . 297622e02 , 1 . 568303e02 , 2 . 302497e 

0026  S03 , - 1 . 540896e02 , 5 . 700617e02 , 1 . 292881e03 , -7 . 922399e02 , 1 . 080333e03 , - 

0027  S  3 . 94 1 08  7  eOl , 2 . 055201 e02 , - 1 . 8531 81 e02 , 3 . 569555e02 , -3 . 656370e01 ,3.01 

0028  S2583e02 , 8 . 903696e01 , -6 . 436587e02 , -2 . 424140e02 , -1 . 041800e03 , 5 . 89817 

0029  S9e02 , 2 . 310479e02 , -5 . 887414e01 , 4 . 001436e01 , -1 . 209943e-02 , 9 . 459898e0 

0030  SO , -1 . 050984e02 , -3 . 745591 e02 , 3 . 563806e02 , -1 . 545264e03 , -6 . B28717e02 , 

0031  Si  681499e02 , 2 . 971388e01 , 6 . 276772eOO, 7 . 309118e01 , -3 . 402882e01 ,3.871 

0032  S370e01 , -1 . 670375e01 , 1 . 915876e03 , 7 . 079673e02 , 1 . 857451e02 , -2 . 732077e 

0033  SOI , -1 . 705171e02 , 5 . 1 1 5862e01 , 1 . 302727e01 , -3 . 776955eOO, -2 . 940332e01 , 

0034  S3 . 510623e-01 , -4 . 633602e02 , 6 . 821298e02 , -2 . 394838e02 , 4 . 549622e02 , -3 . 

0035  $794850e01 , -1 . 617146e02 , 6 . 268821eOO, 1 . 004341e01 , -4 . 002399e00, -4 . 152 

0036  $ 194eOO , 2 . 8031 31 e03 , - 1 . 698787e03 , -4 . 244406e02 , 1 . 998351 e02 , 6 . 192396e 

0037  SOI , -1 . 668931e02 , -9 . 080082e01 , -5 . 963821e-01 , 1 . 524572eOO, -9 . 238670e- 

0038  SOI/ 

0039  data  bm/9 . 933492e04 , 9 . 933492e04 , 3 . 746322e04 , 2 . 457753e04 , 1 . 329481e0 

0040  $4,6. 468820e03 , 3 . 385349e03 , 1 . 61 6258e03 , 7 . 409154e02 , 3 . 641040e02/ 

0041  data  nmax/10/ , berr/O . 0001/ 

0042  c 


BMF  is  declination  of  the  geomagnetic  field 

DIP  is  dip  angle 

B  is  total  field 

BR  is  radial  component 

BP  is  longitudinal  component 

BT  id  latitudinal  component 


0043 

50 

p22=abs (s i n (thet) ) 

0044 

if(p22  .eq.  0.)  p22=l.e-6 

0045 

p21=sqrt(l . -p22*p22) 

0046 

re=6356 . 912*p22*p22* (21 . 3677* . 108*p22*p22) 

0047 

ar=(re+r) /6371 . 2 

0048 

if (thet  . 1 e .  1 ,570796327e0)  go  to  70 

0049 

p21=-p21 

0050 

70 

i f ( j  . eq.  0)  go  to  90 

0051 

ssq=p22*p22 

0052 

ar=ar+(14.288-ssq*(21 . 3677+ . 108*ssq) ) /6371 

0053 

90 

ar=l . /ar 

0054 

c 

n=  2 

0055 

dp22=p21 

0056 

c 

convert  to  east  longitude 

0057 

ph i =ph i j 
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i 

0058 

i f (phi )  92,96,94 

»'  ( 

0059 

92 

ph i =— ph i 

0060 

go  to  96 

0061 

94 

ph i =6 . 2831853e0-ph i 

0062 

96 

sp2=sin(phi) 

0063 

m2=sin(phi) 

0064 

_p2=cos(phi) 

0065 

dp21=-p22 

■  <*# 
vi 

0066 

aor=ar*ar*ar 

0067 

c2=g(2,2) *cp2+g(l ,2) *sp2 

0068 

br=- (aor+aor) * (g (2 , 1) *p21 +c2*p22) 

VV 

0069 

bt=aor* (g(2,l) *dp21+e2*dp22) 

0070 

bp=aor* (g (1,2) *cp2-g (2,2) *sp2) *p22 

0071 

if(nmax  .It.  3)  go  to  260 

f.s* 

0072 

aor=aor*ar 

* . 

0073 

err=berr*sqrt( (bp/p22) **2*br**2+bt*#2) 

0074 

if(bm(3)*aor  .le.  err)  go  to  260 

0075 

sp3=(sp2+sp2) *cp2 

".V 

0076 

cp3=(cp2+sp2) * (cp2-sp2) 

i 

■  V 

0077 

p31=p21*p21-0.333333333e0 

'X*r 

0078 

p32=p21*p22 

0079 

p33=p22*p22 

0080 

dp31=-p32-p32 

•  *  * 

0081 

dp32=p21*p21-p33 

0082 

dp33=-dp31 

0083 

c2=g (3 , 2) *cp2+g (1 , 3) *sp2 

$ 

0084 

c3=g (3 , 3) *cp3+g (2,3)*  sp3 

‘‘V* 

0085 

br=br-3 .0*aor* (g(3, 1) *p31+c2*p32+c3*p33) 

0086 

bt=bt+aor*  (g(3, 1)*dp31«-c2*dp32+c3*dp33) 

0087 

bp=bp-aor * ( (g (3 , 2) *sp2-g (1 , 3) *cp2) *p32+2 . 0* (g (3, 3) *sp3-g (2 , 3) *cp3) 

0088 

$*P33) 

0089 

c 

n=  4 

*\ 

•V ' 

0090 

if(nmax  .It.  4)  go  to  260 

*v> 

0091 

aor=aor*ar 

0092 

if(bm(4)*aor  .le.  err)  go  to  260 

i 

0093 

sp4=sp2*cp3*cp2*sp3 

0094 

cp4=cp2*cp3-sp2*sp3 

0095 

p41=p21*p31-0 . 26666666e0*p21 

0096 

dp41=p21*dp31+dp21*p31-0 . 26666666e0*dp21 

4  ,* 

0097 

p42=p21*p32-0 . 20000000e0*p22 

0098 

dp42=p21*dp32*dp21*p32-0 . 20000000e0*dp22 

0099 

p43=p21*p33 

0100 

dp43=p21*dp33*dp21*p33 

1 

*,>? 
s  v 

0101 

p44=p22*p33 

0102 

dp44=3.0*p43 

0103 

c2=g(4,2) *cp2+g(l ,4) *sp2 

0104 

c3=g(4, 3) *cp3+g (2,4) *sp3 

*  f  ■ 

0105 

c4=g(4,4) *cp4+g(3,4) *sp4 

-T— 

0106 

br=br-4 .0*aor* (g(4 ,1) *p41*c2*p42+c3*p43+c4*p44) 

0107 

bt=bt*aor* (g(4, 1) *dp41*c2*dp42+c3*dp43+c4*dp44) 

»  V 

0108 

bp=bp-aor*((g(4,2)*sp2-g(l,4)*cp2)*p42+2.0*(g(4,3)*sp3-g(2,4)*cp3) 

*>*' 

»v»* 

>.-v 

0109 

$*p43+3.0*(g(4,4) *sp4-g(3,4) *cp4) *p44) 

0110 

if(nmax  .It.  5)  go  to  260 

,*  V 

0111 

aor=aor*ar 

0112 

if(bm(5)*aor  .le.  err)  go  to  260 

0113 

sp5=(sp3*sp3) *cp3 

.  v  *. 

*  Vl 

♦ 

.  1 

0114 

cp5= (cp3+sp3) * (cp3-sp3) 

• 
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0115 

0116 

0117 

0118 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134  c 
0135 
0136 
0137 
0138 
0139 
0140 
0141 
0142 
0143 
0144 
0145 
0146 
0147 
0148 
0149 
0150 
0151 
0152 
0153 
0154 
0155 
0156 
0157 
0158 
0159 
0160 
0161 
0162 
0163 
0164 
0165 
0166 
0167 
0168 
0169 
0170 
0171 


p51=p21 *p41-0 . 25714285e0*p31 
p52=p21 »p42-0 . 22857142e0*p32 
dp51=p21+dp41+dp21*p41-0 . 25714285e0*dp31 
dp52=p21»dp42+dp21*p42-0 . 228571 42e0*dp32 
p53=p21*p43-0 . 14285714e0*p33 
dp53=p21*dp43+dp21*p43-0 . 14285714e0*dp33 
p54=p21*p44 

dp54=p21*dp44+dp21*p44 

p55=p22*p44 

dp55=4 . 0*p54 

c2=g(5, 2) *cp2+g (1 , 5) *sp2 
c3=g (5 , 3) *cp3+g (2 , 5) *sp3 
c4=g(5,4) *cp4+g(3, 5) *sp4 
c5=g(5 , 5) *cp5+g(4,5) *sp5 

br=br-5 .0*aor* (g (5,1) *p51+c2*p52+c3*p53+c4*p54+c5*p55) 
bt=bt+aor* (g(5, 1) *dp51+c2*dp52+c3*dp53+c4*dp54+c5*dp55) 
bp=bp-aor* ( (g(5,2) *sp2-g(l ,5) *cp2) *p52+2.0* (g(5,3) *sp3-g(2,5)*cp3) 
$*p53+3 .0* (g(5,4) *sp4-g(3,5) *cp4) *p54+4 .0* (g(5,5) *sp5-g(4,5)*cp5) *p 
$55) 
n=  6 

if(nmax  .It.  6)  go  to  260 
aor=aor*ar 

if(bm(6)*aor  .le.  err)  go  to  260 
sp6=sp2*cp5+cp2*sp5 
cp6=cp2*cp5-sp2*sp5 
p61=p21*p51-0 . 25396825e0*p41 
dp61=p21*dp51+dp21*p51-0.25396825e0*dp41 
p62=p21*p52-0 . 23809523e0*p42 
dp62=p21*dp52+dp21*p52-0 . 23809523e0*dp42 
p63=p21*p53-0 . 19047619e0*p43 
dp63=p21*dp53+dp21*p53-0 . 19047619e0*dp43 
p64=p21*p54-0. Illllllle0*p44 
dp64=p21 *dp54+dp21*p54-0. 1111111 Ie0*dp44 
p65=p21*p55 

dp65=p21*dp55+dp21*p55 

p66=p22*p55 

dp66=5 . 0*p65 

c2=g(6,2)*cp2+g(l,6)*sp2 
c3=g(6,3) »cp3+g (2 , 6) *sp3 
c4=g(6,4) *cp4+g(3, 6) *sp4 
c5=g(6,5) *cp5+g(4,6) *sp5 
c6=g(6,6) *cp6+g(5,6) *sp6 

br=br-6 . 0*aor* (g (6, 1) *p61+c2*p62+c3*p63+c4*p64+c5*p65+c6*p66) 
bt=bt+aor* (g (6, 1) *dp61+c2*dp62+c3*dp63+c4*dp64+c5*dp65+c6*dp66) 
bp=bp-aor* ( (g (6 , 2) *sp2-g(l , 6) *cp2) *p62+2 .0* (g(6,3) *sp3-g(2,6) *cp3) 
$*p63+3 .0* (g (6,4) *sp4-g (3,6) *cp4) *p64+4 .0* (g(6, 5) *sp5-g(4,6) *cp5) *p 

$65+5 . 0* (g (6 , 6) *sp6-g (5 , 6) *cp6) *p66) 
if(nmax  .It.  7)  go  to  260 
aor=aor*ar 

i f  (bm(7) «aor  .le.  err)  go  to  260 

sp7=(sp4+sp4) *cp4 

cp7= (cp4+sp4) * (cp4-sp4) 

p71=p21 *p61-0 . 25252525e0*p51 

dp71=p21 »dp61+dp21*p61-0. 25252525e0*dp51 

p72=p21 *p62-0 . 24242424e0*p52 

dp72=p21*dp62+dp21*p62-0 . 24242424e0»dp52 

p73=p21*p63-0 . 21212121e0*p53 
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0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 

0181 

0182 

0183 
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0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195  c 
0196 
0197 
0198 
0199 
0200 
0201 
0202 
0203 
0204 
0205 
0206 
0207 
0208 
0209 
0210 
0211 
0212 
0213 
0214 
0215 
0216 
0217 
0218 
0219 
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0221 
0222 
0223 
0224 
0225 
0226 
0227 
0228 


dp73=p21*dp63*dp21*p63-0.21212121e0*dp53 
p74=p21*p64-0. 16161616e0*p54 
dp74=p21*dp64+dp21*p64-0. 16161616e0*dp54 
p75=p21*p65-0  09090909e0*p55 
dp75=p21 *dp65*dp21 *p65-0 . 09090909e0*dp55 
p76=p21*p66 

dp76=p21*dp66*dp21*p66 

p77=p22*p66 

dp77=6.0*p76 

c2=g(7 ,2) *cp2+g(l ,7) *sp2 
c3=g(7,3) *cp3-*-g(2,7) *sp3 
c4=g(7,4)*cp4*g(3J7)*sp4 
c5=g (7 , 5) *cp5*g (4 , 7) *sp5 
c6=g (7 , 6) *cp6*g (5 , 7) *sp6 
c7=g(7,7)*cp7+g(6,7)*sp7 

br=br-7 .0*aor* (g(7 , 1) *p71*c2*p72+c3*p73+c4*p74+c5*p75+c6*p76*c7*p7 
$7) 

bt=bt+aor* (g(7 , 1) *dp71+c2*dp72+c3*dp73«-c4*dp74*c5*dp75+c6*dp76+c7* 
$dp77) 

bp=bp-aor* ( (g(7 ,2) *sp2-g(l ,7) *cp2) *p72+2 .0* (g(7,3) *sp3-g(2, 7) *cp3) 
S»p73+3.0*(g(7,4)*sp4-g(3,7)*cp4)*p74+4.0*(g(7,5)*sp5-g(4,7)*cp5)*p 
$75+5 .0* (g (7 , 6) *sp6-g (5,7) *cp6) *p76+6 .0* (g (7 ,7) *sp7-g(6, 7) *cp7) *p77 

S) 

n=  8 

if(nmax  It.  8)  go  to  260 
aor=aor*ar 

if(bm(8)*aor  .le.  err)  go  to  260 
sp8=sp2*cp7*cp2*sp7 
cp8=cp2*cp7-sp2*sp7 
p81=p21 *P71 -0 . 25174825eO*p61 
dp81=p21*dp71*dp21*p71-0. 25174825e0*dp61 
p82=p21*p72-0 . 24475524e0*p62 
dp82=p21*dp72+dp21*p72-0 . 24475524e0*dp62 
p83=p21.p73-0.22377622e0*p63 
dp83=p21*dp73*dp21*p73-0.22377622eO*dp63 
p84=p21*p74-0. 18881118e0*p64 
dp84=p21*dp74+dp21*p74-0. 18881 11 8e0*dp64 
p85=p21*p75-0. 13986013e0*p65 
dp85=p21*dp75+dp21*p75-0 . 13986013e0*dp65 
p86=p21*p76-0 . 07692307e0*p66 
dp86=p21*dp76+dp21*p76-0 . 07692307e0*dp66 
p87=p21*p77 

dp87=p21»dp77*dp21«p77 

p88=p22*p77 

dp88=7.0*p87 

c2=g(8,2)*cp2*g(l,8)*sp2 
c3=g(8,3) *cp3*g(2, 8) *sp3 
c4=g (8 , 4) *cp4+g (3 , 8) *sp4 
c5=g (8 , 5) *cp5+g (4,8) *sp5 
c6=g (8,6) *cp6+g(5, 8) *sp6 
c7=g(8,7)*cp7*g(6,8)*sp7 
c8=g(8 , 8) *cp8*g(7 , 8) *sp8 

br=br-8  0*aor* (g(8, 1) *p81*c2*p82*c3*p83+c4*p84+c5*p85*c6*p86+c7*p8 
17+c8*p88) 

bt=bt+aor* (g(8, 1) *dp8l ♦c2*dp82*c3*dp83+c4*dp84»c5*dp85+c6*dp86*c7* 
Jdp87*c8*dp88) 

bp=bp-aor*((g(8,2)*sp2-g(l,8)*cp2)*p82'*2.0*(g(8,3)*sp3-g(2,8)*cp3) 
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0231 
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$*p83+3.0*(g(8,4)*sp4-g(3,8)*cp4)*p84+4.0* (g(8,5)*sp5-g(4,8)*cp5)*p 
$85+5.0* (g(8,6) *sp6-g(5,8) *cp6) *p86+6 .0* (g(8,7) *sp7-g(6,8) *cp7)*p87 
$+7.0* (g(8,8) *sp8-g(7,8) *cp8) *p88) 
if(nmax  .It.  9)  go  to  260 
aor=aor *ar 

i  f  (bin (9)  *aor  .le.  err)  go  to  260 
sp9=(sp5+sp5) *cp5 
cp9=(cp5+sp5) * (cp5-sp5) 
p91=p21*p81-0 . 25128205e0*p71 
dp91=p21*dp81+dp21*p81-0 . 25128205e0*dp71 
p92=p21*p82-0.24615384e0*p72 
dp92=p21*dp82+dp21*p82-0.24615384e0*dp72 
p93=p21*p83-0 . 23076923e0*p73 
dp93=p21 *dp83+dp21 *p83-0 . 23076923e0*dp73 
p94=p21*p84-0 . 20512820e0*p74 
dp94=p21*dp84+dp21*p84-0 . 20512820e0*dp74 
p95=p21*p85-0 . 16923076e0*p75 
dp95=p21*dp85+Hp21*p85-0 . 16923076e0*dp75 
p96=p21*p86-0 . 12307692e0*p76 
dp96=p21*dp86+dp21*p86-0 . 12307692e0*dp76 
p97=p21 *p87-0 . 06666666e0*p77 
dp97=p21*dp87+dp21*p87-0 ,06666666eO*dp77 
p98=p21*p88 

dp98=p21*dp88+dp21*p88 

p99=p22*p88 

dp99=8 . 0*p98 

c2=g(9, 2) *cp2+g (1 , 9) *sp2 
c3=g(9,3) *cp3+g(2,9)*sp3 
c4=g (9 , 4) *cp4+g (3 , 9) *sp4 
c5=g (9 , 5) *cp5+g (4 , 9) *sp5 
c6=g(9,6) *cp6+g(5,9)*sp6 
c7=g(9, 7) *cp7+g(6,9)*sp7 
c8=g(9,8) *cp8+g(7,9) *sp8 
c9=g (9 , 9) *cp9+g (8 , 9) *sp9 

br=br-9  0*aor* (g (9, 1) *p91+c2*p92+c3*p93+c4*p94+c5*p95+c6*p96+c7*p9 
$7+c8*p98+c9*p99) 

bt=bt+aor* (g(9, 1) +dp91+c2*dp92+c3*dp93+c4*dp94+c5*dp95+c6*dp96+c7* 
$dp97+c8*dp98+c9*dp99) 

bp=bp-aor * ((g (9 , 2) *sp2-g (1,9) *cp2) *p92+2 .0* (g (9,3) *sp3-g(2,9) *cp3) 
S*p93+3  0* (g(9,4) *sp4-g(3,9) *cp4) *p94+4 .0* (g(9,5) *sp5-g(4,9) *cp5) *p 
$95+5.0* (g (9 , 6) *sp6-g (5,9) *cp6) *p96+6 .0* (g(9,7) *sp7-g(6,9) *cp7)*p97 
$+7.0* (g(9,8) *sp8-g(7,9) *cp8) *p98+8 . 0* (g(9, 9) *sp9-g(8,9) *cp9) *p99) 
n=10 

if(nmax  .It.  10)  go  to  260 
aor=aor»ar 

i f  (bm(10) *aor  .le.  err)  go  to  260 
spl0=sp2*cp9+cp2*sp9 
cpl0=cp2*cp9-sp2*sp9 
pl01=p21*p91-0 . 25098039e0»p81 
dpl01=p21*dp91+dp21*p91-0 . 25098039e0*dp81 
pl02=p21 *p92-0 . 24705882eO*p82 
dpl02=p21*dp92+dp21*p92-0 . 24705882e0*dp82 
pl03=p21 *p93-0 . 23529411e0*p83 
dpl03=p21*dp93+dp21*p93-0. 2352941 Ie0*dp83 
pl04=p21*p94-0 . 21568627eO*p84 
dpl04=p21*dp94+dp21*p94-0.21568627e0*dp84 
pl05=p2l*p95-0.18823529e0*p85 
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dpl05=p21-*dp95+dp21*p95-0 . 18823529e0*dp85 

pl06=p21*p96-0. 15294117e0*p86 

dpl06=p21*dp96+dp21*p96-0.15294117e0*dp86 

pl07=p21*p97-0 . 10980392e0*p87 

dpl07=p21*dp97*dp21*p97-0 . 10980392e0*dp87 

pl08=p21*p98-0.05882352e0*p88 

dpl08=p21*dp98*dp21*p98-0. 05882352e0*dp88 

pl09=p21*p99 

dpl09=p21*dp99+dp21*p99 

plOlO=p22*p99 

dpl010=9.0*pl09 

c2=g(10,2) *cp2+g(l, 10) *sp2 

c3=g(10,3) *cp3+g(2,10)*sp3 

c4=g (10, 4) *cp4+g(3, 10) *sp4 

c5=g(10, 5) *cp5+g(4, 10) *sp5 

c6=g (10,6) *cp6+g(5,10) *sp6 

c7=g(10,7) *cp7+g(6,10) *sp7 

c8=g(10,8) *cp8+g(7,10) *sp8 

c9=g (10, 9) *cp9+g (8 , 10) *sp9 

cl0=g(10, 10) *cpl0+g(9, 10) *splO 

br=br-10 .0*aor* (g(10, 1) *pl01+c2*pl02+c3*pl03+c4*pl04*c5*pl05+c6*pl 
$06*c7*pl07+c8*pl08+c9*pl09+cl0*pl010) 

bt=bt+aor* (g(10, 1) *dpl01+c2*dpl02+c3*dpl03*c4*dpl04+c5*dpl05+c6*dp 
1106+c7*dpl07*c8*dpl08+c9*dpl09+cl0*dpl010) 
bp=bp-aor* ( (g (10, 2) *sp2-g(l , 10) *cp2) *pl02*2 .0* (g(10, 3) *sp3-g (2,10) 
I*cp3) *pl03+3 . 0* (g (10, 4) *sp4-g (3 , 10) *cp4) *pl04+4 . 0* (g (10, 5) *sp5-g (4 
S , 10) *cp5) *pl05*5  0* (g (10,6) *sp6-g (5, 10) *cp6) *pl06*6 . 0* (g (10, 7) *sp7 
J-g(6, 10)  *cp7)  *pl07*7 .0*  (g(10,8)  *sp8-g(7, 10)  *cp8)  *pl08-»8.0*(g(10,9) 
S*sp9-g (8 , 10) *cp9) *pl09*9 .0* (g (10, 10) ♦splO-g (9,10) ♦cplO) *pl010) 

260  bp=bp/p22*l . e-5 
bt=bt*l . e-5 
br=br*l .e-5 

b=sqrt(br*br*bt*bt+bp*bp) 

bh=sqrt (bt*bt+bp*bp) 

bmf =3 . 141592654eO-acos  (bt/bh) 

if (bp  .It.  0.)  bmf=-bmf 

dip=acos(bh/b) 

if(br  .gt.  0.)  dip=-dip 


0001 

subroutine  prof i n(l u, type,maxhts, npr i nt , nrspec, Imax, h 1 i st,a logen) 

0002 

c 

0003 

c 

Reads  ionospheric  profiles 

0004 

c 

type=l :  electron  and  ion  densities 

0005 

c 

2:  collision  frequencies 

0006 

c 

0007 

i nteger  type 

0008 

character*80  bed 

0009 

d i mens  ion  h 1 i st  (maxhts) , a  1 ogen  (maxhts, 3) , en (3) 

0010 

if (type  .ne.  2)  then 

0011 

read (lu, 1010)  bed 

0012 

if(nprint  .gt.  1)  print  1011, bed 

0013 

end  i  f 

0014 

do  202  l=l,maxhts+l 

0015 

read ( 1 u , 1020, end=900)  ht,en 

0016 

if(ht  .It.  0.)  then 

0017 

lmax=l-l 

0018 

return 

0019 

end  i  f 

0020 

i f ( 1  .ne.  1  .and.  ht  .ge.  hlist(l-l))  then 

0021 

print  *, ’ERROR  PROFIN:  Profile  heights  out  of  order* 

0022 

go  to  999 

0023 

end  i  f 

0024 

hi ist(l)=ht 

0025 

if (type  .eq.  1  .and.  nrspec  .eq.  3)  en(3)=en(2)-en(l) 

0026 

if(nprint  .gt.  1)  print  1021 , ht, (en (k) , k=l , nrspec) 

0027 

do  201  k=l, nrspec 

0028 

201 

alogen(l , k)=a 1 og(amaxl (en (k) ,1 .e-20)) 

0029 

202 

continue 

0030 

print  *, ’ERROR  PROF IN:  Too  many  heights  in  profile* 

0031 

go  to  999 

0032 

900 

print  *, ’ERROR  PROFIN:  Profile  input  not  properly  terminated’ 

0033 

999 

lmax=-l 

0034 

return 

0035 

1010 

f ormat(a80) 

0036 

1011 

format(/lx,a80) 

0037 

1020 

format(f7.2,4x,3el0.2) 

0038 

1021 

format (f 8 . 2, 4x, Ip3el0 .2) 

0039 

end 
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subroutine  qartic (q,b3, b2,bl,b0, debug, newq) 
implicit  real  *8  (a-h,o-z) 

comp  I  ex*  16  b3, b2, bl , bO, q,b3sq, h,  i ,g, hpr  ime,gpr ime,sqroot, 

S  pi , p21cbert0,cbertl , cbert2,omegal , omega2, 

S  rootp, rootq, rootr, fnc ton, c temp, dfdq,dq 

integer  debug 
dimension  diff(4),q(4) 
data  omegal/ (-5 . d-1 ,  8 . 660254038d-l) / 
data  omega2/ (-5 . d-1 , -8 . 660254038d-l) / 
data  tol/1 .d-06/, imax/5/ 
c 

iagain=0 

if (newq  .eq.  1)  go  to  30 
newq=l 

10  b3sq=b3**2 

h=b2-b3sq 

i =b0- (4 . d0,0 . dO) *b3*bl+ (3 . d0,0 . dO) *b2**2 
g=bl+b3* ( (-3 . dO , 0 . dO) *b2+ (2 . dO ,0 . dO) *b3sq) 
hpr ime=-i /(12.d0,0.d0) 

gpr ime=-g**2/ (4 .d0,0 .  dO) -h* (h**2+ (3 .dO.O.dO) *hpr ime) 
c 

sqroot=cdsqrt(gpr ime* *2+ (4 .d0,0.d0) *hpr ime* *3) 
pl=(- . 5d0, 0 . dO) * (gpr ime-sqroot) 
p2=(- .5d0,0 .dO) * (gpr ime* sq root) 
if(cdabs(pl)  .It.  cdabs(p2))  pl=p2 
cbertO=cdexp(cdl og(pl)/(3 .dO^O  dO)) 
cbertl=omegal*cbertO 
cbert2=omega2*cbert0 
c 

rootp=cdsqrt(cbertO-hprime/cbertO-h) 
rootq=cdsqrt(cbertl-hpr ime/cbertl-h) 
rootr=cdsqrt(cbert2-hpr ime/cbert2-h) 
if(cdabs(g)  .gt.  l.d-30)  then 
s i gn=- rootp* rootq*rootr* (2 . d0,0 . dO) /g 
if (sign  .It.  O.dO)  rootr=-rootr 
end  i  f 

q(l)=+ rootp* rootq+rootr-b3 
q(2)=+rootp-rootq-rootr-b3 
q(3)=-rootp*rootq-rootr-b3 
q(4) =-rootp-rootq+rootr-b3 
c 

30  if(debug  .gt.  2)  print  100,b3, b2,bl,b0 
do  60  n=l , 4 
do  40  i ter=l , imax 

fncton=(((q(n)+(4.d0,0.d0)*b3)*q(n)+(6.d0,0.d0)*b2)*q(n) 

$  +(4.60,0. dO) *bl)*q(n)*bO 

dfdq= ( ( (4 . dO , 0 . dO) *q (n) + (12 . dO , 0 . dO) *b3) *q (n) 

$  + (12 . dO , 0 . dO) *b2) *q (n) + (4 . dO, 0 . dO) *bl 

dq=-f ncton/df dq 
q(n)=q(n)+dq 
testdq=cdabs (dq/q (n) ) 
if(testdq  .le.  tol)  go  to  60 
40  continue 

if (i again  .eq.  1)  then 

f ncton=( ((q(n) ♦ (4 .d0,0 .dO) *b3) *q(n) ♦ (6 ,d0,0 .dO) *b2) *q(n) 

S  +(4.d0,0.d0)*bl)*q(n)+b0 

print  101,n,q(n) ,fncton,dq, i again 


i aga i n=l 
go  to  10 
end  i  f 
conti nue 

1=0 

do  80  m=2,4 
do  80  n=m,4 

i f (d i mag (q (n) )  . gt .  0 . dO)  go  to  80 
1  =  1+1 

ctemp=q(n) 
q (n) =q (m-1) 
q(m-l)=ctemp 
conti nue 

i f ( I  . eq .  2)  go  to  99 
do  81  n=l , 4 

a  ngq=cdang (q (n) ) *57 . 29577951 3d0 
if(angq  .It.  135, dO)  angq=angq+360.d0 
di ff (n)=dabs(angq-315.d0) 
do  82  nm=2,4 
do  82  n=nm,4 

if(diff(n)  .gt.  diff(nm-l))  go  to  82 

temp=d i f f (n) 

d i ff (n)=di f f (nm-1) 

diff (nm-l)=temp 

ctemp=q(n) 

q(n)=q(nm-l) 

q (nm-1) =ctemp 

conti nue 

return 

format(/’  In  QARTIC:  b”s  =’ ,4(lpel3.4,el2.4)) 
f ormat(8h  q  root  ,il,2h  =,lp2el3.5,3x,10hfunction  =,2el3.5,3x 
4hdq  =,2el3 . 5,3x, 8h i aga i n  =,il) 


ft 


0001 

subroutine  rbars 

0002 

implicit  real  *8  (a-h,o-z) 

0003 

c 

0004 

include  ’commonl . for ’ 

0020 

include  ’common2 . for  * 

0034 

include  ’common3 . for  * 

0059 

c 

0060 

complex*16  ngsq, sqroot, ratio, ikc,exd,exdsq,zl,z2l 

0061 

$  pO, hlO, h20, hlprmO, h2prm0, caphl0,caph20, 

0062 

$  pd, hid, h2d, hlprmd, h2prmd,caphld,caph2d, 

0063 

$  alst,a2nd,a3rd , a4th, al,a2,a3,a4,fl,f2 

0064 

real*8  kvraot, kvratt, ndsq , nOsq 

0065 

c 

0066 

ngsq=dcmp 1 x (db 1 e (epsr) , -db 1 e (s i gma) / (omega*8 . 85434d-12) ) 

0067 

sqroot=cdsqrt(ngsq-ssq) 

0068 

c 

0069 

i f (d i mag (theta)  .It.  -10. dO  .or.  alpha  .  eq.  0.)  go  to  20 

0070 

if(d  .eq.  0.)  go  to  10 

0071 

c 

0072 

kvraot=dexp  (d 1 og  (wn/a 1 pha) /3 . dO) 

0073 

kvratt=kvraot**2 

0074 

avrkot=l .dO/kvraot 

0075 

a vrktt=avrkot**2*0 . 5d0 

0076 

nOsq=l . -a  1 pha*h 

0077 

rati o=nOsq/ngsq*sqroot 

0078 

pO=kvratt* (nOsq-ssq) 

0079 

call  mdhnk  1  (pO, hlO, h20, hlprmO, h?;prmO, theta rb  1’) 

0080 

caphlO=hlprmO+avrktt*hlO 

0081 

caph20=h2prm0*avrktt*h20 

0082 

alst=caph20-zmp 1 x i *rati o*kvraot*h20 

0083 

a2nd=caphl0-zmpl xi *rati o*kvraot*hl0 

0084 

a3rd=h2prmO-zmp 1 x i *k vraot*sqroot*h20 

0085 

a4th=hlprmO-zmp 1 x i *kvraot*sqroot*hlO 

0086 

ndsq=l . -a  1 pha* (h-d) 

0087 

pd=kvratt* (ndsq-ssq) 

0088 

call  mdhnk 1  (pd, hid, h2d, hlprmd , h2prmd, theta rb  2’) 

0089 

caphld=hlprmd*avrktt*hld 

0090 

caph2d=h2prmd+avrktt*h2d 

0091 

f I=h2d*a2nd-hld*alst 

0092 

f 2=h2d*a4th-hld*a3rd 

0093 

al=c*ndsq*f 1 

0094 

a2=zmp 1 x i *avrkot* (caphld*alst-caph2d*a2nd) 

0095 

a3=zmp 1 x i *avrkot* (h2prmd*a4th-hlprmd*a3rd) 

0096 

a4=c*f2 

0097 

rbarll=(al-a2)/(al+a2) 

0098 

rbar22= (a3+a4) / (a4-a3) 

0099 

hg=exp(- ,5*a lpha*d) *(h20*a2nd-hl0*alst)/f 1 

0100 

normll=f l*f 1 

0101 

norm22=f 2*f 2 

0102 

norml2=f l*f 2 

0103 

return 

0104 

c 

0105 

10 

rbarll=(ngsq*c-sqroot) / (ngsq*c+sqroot) 

0106 

rbar22=(c-sqroot) / (c+sqroot) 

0107 

hg=zone 

0108 

normll=(-2 . 124292958d0, 0. dO) 

0109 

norm22=normll 

RBARS 

0110 

norml2=normll 

0111 

return 

0112 

c 

0113 

c 

flat  earth 

0114 

20 

i kc=dcmp 1 x (0 . dO , -wn) *c 

0115 

exd=cdexp(i kc*d) 

0116 

exdsq=exd*exd 

0117 

zl=(ngsq*c-sqroot) / (ngsq*c+sqroot) 

0118 

z2=(c-sqroot) / (c+sqroot) 

0119 

rbarll=zl*exdsq 

0120 

rbar22=z2*exdsq 

0121 

hg=exd* (zone+zl) /(zone+rbarll) 

0122 

normll=(zone+rbarll) * (zone+rbarll) /exdsq 

0123 

norm22=(zone+rbar22) * (zone+rbar22) /exdsq 

0124 

norml2= (zone+rbarll) *  (zone+rbar22) /exdsq 

0125 

return 

0126 

end 

0001 

subrout i ne  rec vr (t 1 ng, tc 1 t, xtr ,  rho, r 1 ng, rc 1 t) 

0002 

c 

0003 

c 

Returns  coordinates  of  a  point  which  is  at  a  specified  great 

0004 

c 

circle  distance  and  bearing  angle  from  the  input  point 

0005 

c 

0006 

c 

Input:  TLNG  is  longitude  of  transmitter 

0007 

c 

TCLT  is  co-latitude  of  transmitter 

0008 

c 

XTR  is  geographic  bearing  angle  of  receiver 

0009 

c 

RHO  is  great  circle  distance  to  the  receiver 

0010 

c 

0011 

c 

Output:  RLNG  is  longitude  of  receiver 

0012 

c 

RCLT  is  co- latitude  of  receiver 

0013 

c 

0014 

c 

All  coordinates,  RHO  and  XTR  are  in  radians 

0015 

c 

Sign  convention  is  ♦  for  West  and  North 

0016 

c 

0017 

data  p i /3 . 14159265eO/ , twop i /6 . 28318531e0/ 

0018 

c 

0019 

red uce(arg)=sign (ami nl(abs(arg) ,1 .) ,arg) 

0020 

c 

0021 

etc  1 t=cos (tc 1 1) 

0022 

stc 1 t=si n (tc 1 t) 

0023 

br=xtr 

0024 

gcd=rho 

0025 

c 

0026 

if(abs(br)  .It.  twopi)  go  to  2 

0027 

br=amod (br , twop i ) 

0028 

2 

if(br  .ge.  0.)  go  to  3 

0029 

br=br*twop i 

0030 

3 

if (ged  .It.  pi)  go  to  5 

0031 

gcd=twopi -ged 

0032 

br=br+p i 

0033 

if(br  .ge.  twopi)  br=br-twopi 

0034 

5 

if(br  ,le.  l.e-6)  go  to  10 

0035 

i f  (abs  (br-p i )  ,le.  l.e-6)  go  to  14 

0036 

i f  (abs  (ged-p i )  .  le.  l.e-6)  go  to  14 

0037 

cgcd=cos (ged) 

0038 

sgcd=s i n  (ged) 

0039 

crc 1 t=ctc 1 t*cgcd+stc 1 t*sgcd*cos (br) 

0040 

sre 1 t=sqrt(l . -crc 1 t**2) 

0041 

rc 1 t=acos ( reduce (crc 1 t) ) 

0042 

de 1 ta=acos (reduce ( (cgcd-ctc 1 t*crc 1 t) / (stc 1 t*src 1 t) ) ) 

0043 

i f (br  .It.  pi)  delta=-delta 

0044 

r 1 ng=t 1 ng+de 1 ta 

0045 

go  to  20 

0046 

c 

0047 

c 

receiver  is  due  north,  south  or  on  opposite  longitude 

0048 

10 

rc 1 t=tc 1 t-ged 

0049 

i f (rc 1 t  .It.  0 . )  go  to  12 

0050 

11 

r  1  ng=t !  <g 

0051 

crc 1 t=cos(rc 1 t) 

0052 

sre 1 t=s i n (rc 1 t) 

0053 

go  to  99 

0054 

12 

rc 1 1=- rc  1 1 

0055 

13 

r 1 ng=tl ng+pi 

0056 

crc 1 t=cos(rc 1 t) 

0057 

sre 1 t=s i n(rc 1 t) 
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0058 

go  to  20 

0059 

14 

rclt=tclt+gcd 

0060 

if (rclt  .It.  pi) 

go 

to 

11 

0061 

re  1 t=twopi -rclt 

0062 

go  to  13 

0063 

c 

0064 

20 

i f (r 1 ng  ,gt.  pi) 

go 

to 

21 

0065 

if (ring  .It.  -pi) 

go 

to 

22 

0066 

go  to  99 

0067 

21 

r 1 ng=r 1 ng-twop i 

0068 

go  to  99 

0069 

22 

r lng=rlng+twopi 

0070 

c 

0071 

99 

return 

0072 

end 
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subroutine  rplyn* 
implicit  real  *8  (a-h,o-*) 
c 

include  ’co*mon2 . f or ' 
include  'coeaion3 .  for  * 
c 

comp  I  ex* 16  lgmtrx(30,4) , prod, 1 1 istl,tl ist2 
complex*  8  stheta 
real *4  dst(30) 
integer  use (30) 
c 

m*l 

10  if(m  le  30  and  tlist(l,m)  gt  0  )  then 
thetarstl i st(l ,m) 
theta i *tl i st(2,m) 
c =c dc os (theta* idtr) 
esqae  *c 

s=cds  n(theta*idtr) 
ssqas*s 
call  nteg 
do  12  n= 1 ,4 

12  l gmtr * (m, n) = I ogr s (n) 

ad  j  f I g*l 
m=»*  1 
go  to  10 
end  i  f 
max=m- 1 

if(imax  le  1)  then 

pr  i  nt  *,’EJ»08  RPIYI8I  Insufficient  tlist’ 
stop 

else 

jmaiiamirtOf  ,  nr  t  I  it ) 

ad j  # i gaO 
I-Sto-  1 

end  f 

c 

entry  .  >pc  y 

c  u  stance  * '  qm  tneta  to  t  st  ang 
strata*  \*eta 
dc  ?4  =  1  *a> 

j«* '  1  ,  *  1 

?4  dst  l,*sqrt''  rea 1  fstf*etay  - 1  st(l  lj;**2* 

I  a  mag(sthetaj - t i  at!?  li/**2, 

c  Orde'  t  st  ang  ei  accord  ng  to  distance 

ceM  so'tr  (dst ,  me>  ,  use  me*  .  1  max; 
c  Use  only  nrt 1st  angles 
do  50  n«l  4 

ogr  s'  r  *0 

dc  46  ,1  =  1  j  me  • 

1  -u se ( .  ; , 

1 1 ' stl *dcmp i ■ ( dto i e( 1 1 ■ st O  1|  dto’efti  *112,  ' , 
pr  od* /one 
dc  44  j  2»1 , jme* 

2  =  u*e  '  j? , 

4  i  l  2;  tie* 

t  s  1 2  «dr  a*  >  It  *  t '  1  /  dt  e  |  *  » 1  1  2 

pr odapr od* ' t'eta - t  St?) / ( t  Stl -t  st?; 
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0095 

end  if 

0096 

44 

continue 

0097 

45 

logrs(n)  =  logrs(n) 4-prod*  lgmtrx(il,n) 

0098 

50 

rs(n)=cdexp( logrs(n)) 

0099 

return 

0100 

end 
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subroutine  savemc 

This  routine  writes  the  mode  parameters  out  to  the  logical  unit 
def i ned  by  LUNIT7 . 

include  'commonl .for ’ 
include  ’common2.for ’ 

wr i te ( I un i t7 , 100)  r ho , f r eq , az i m , cod i p , magf I d , s i gma , epsr , hprout 
do  10  m=l, modes 

write(lunit7,101)  tp(m) ,nterm(m) ,t  term(l,m),t  term(2,m), 

S  tp(m) ,nterm(m) ,t  term(3,m),t  term(4,m) 

write(lunit7,102) 
return 

format(’ r * ,f 7 .3, *  f\f8.4,’  a\f8.3,’  c’.fe.S,’  m’,el0.3, 

S  *  s’ ,lpel0.3,  *  e ’  ,0pf 5 . 1 ,  ’  t\f5.1) 

formate  1 ’ (0p2f 9 . 5 , i 1 , Ip4el5 . 8/ ’2 ’ ,0p2f 9 .5, il,lp4el5.8) 
formate  ’) 
end 


0001 

subroutine  sortr (array, nra, index, nr i ,  i  i ,  j j) 

0002 

c 

0003 

c 

algorithm  347, r .c . si ng 1 eton,commun i cati ons  of  the  acm, vl2, n3,mar69 

0004 

c 

sorts  array  into  order  of  increasing  value,  from  index  ii  to  jj 

0005 

c 

also  orders  index  simultaneously  if  nri  gt  1 

0006 

c 

the  only  arithmetic  operatic  on  array  is  subtraction 

0007 

c 

the  user  should  consider  the  possibility  of  integer  overflow 

0008 

c 

arrays  i u (k)  cud  i 1 (k)  permit  sorting  up  to  2**(k*l)-l  elements 

0009 

c 

0010 

dimension  array (1) , index(l) , iu(36) , i 1 (36) 

0011 

i f  ( j  j  gt.  nra)  print  *,’warning  from  sortr:  jj  gt  nra* 

0012 

m=l 

0013 

i  =  i  i 

0014 

j=jj 

0015 

5 

if(i  . ge.  j)  go  to  70 

0016 

10 

k=i 

0017 

ij=(i*j)/2 

0018 

t=array (i j) 

0019 

if (nri  . le.  1)  go  to  15 

0020 

n= i ndex ( i j) 

0021 

15 

if  (array  (i)  le.  t)  go  to  20 

0022 

array ( i j )=array  ( i ) 

0023 

array(i)=t 

0024 

t=array ( i j ) 

0025 

i f (nr i  . 1 e .  1)  go  to  20 

0026 

i ndex ( i j )=i ndex ( i ) 

0027 

i ndex ( i )=n 

0028 

n=i ndex(i j) 

0029 

20 

l=j 

0030 

if(array(j)  .ge.  t)  go  to  40 

0031 

array (ij)*array(j) 

0032 

array(j)=t 

0033 

t=array (i j) 

0034 

i f (nr i  . le.  1)  go  to  25 

0035 

i ndex ( i j )=i ndex (j ) 

0036 

i ndex (j ) =n 

0037 

n= i ndex ( i j ) 

0038 

25 

if (array (i)  le.  t)  go  to  40 

0039 

array ( i j )=array ( i ) 

0040 

array ( i )=t 

0041 

t=array ( i j ) 

0042 

i f  (nr  i  . 1 e .  1)  go  to  40 

0043 

i ndex ( i j ) = i ndex ( i ) 

0044 

i ndex ( i ) =n 

0045 

n= i ndex ( i j ) 

0046 

go  to  40 

0047 

30 

array ( 1 )=array  (k) 

0048 

array (k) =tt 

0049 

i  f  (nr  i  le  1)  go  to  40 

0050 

t ndex ( 1 )  =  i ndex  (k) 

0051 

i ndex (k) *nn 

0052 

40 

1  =  1-1 

0053 

i  f  (array  (1)  gt  t)  go  to  40 

0054 

tt=array  ( 1 ) 

0055 

i  f  (nr  i  le.  1)  go  to  50 

0056 

nn= i ndex  ( 1 ) 

0057 

50 

k=k*l 
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0058 

if(array(k)  .It.  t)  go  to 

50 

0059 

i  f  (k  .  1  e .  1 )  go  to  30 

0060 

if(l-i  .!e.  j-k)  go  to  60 

0061 

i  1  (m)  =  i 

0062 

i  u(m)  =  l 

0063 

i=k 

0064 

♦  1 

0065 

go  to  80 

0066 

60 

i  1  (m) =k 

0067 

i  u  (m)  =  j 

0068 

j  =  l 

0069 

di-in*  1 

0070 

go  to  80 

0071 

70 

m=m- 1 

0072 

i f (m  .  eq.  0)  return 

0073 

i=i 1 (m) 

0074 

j  =  i  u  (m) 

0075 

80 

i  f  (j -  i  ge  11)  go  to  10 

0076 

i  f  ( i  eq  i  i )  go  to  5 

0077 

i  =  i-l 

0078 

90 

i  =  i  + 1 

0079 

i f ( i  .  eq .  j )  go  to  70 

0080 

t=array ( i *1) 

0081 

if  (nr i  .  le.  1)  go  to  95 

0082 

n=i nde*(i *1) 

0083 

95 

i f  (array  ( i )  le.  t)  go  to 

90 

0084 

k=  i 

0085 

100 

array  (k*l)«ar ray (k) 

0086 

i f (nr i  le.  1)  go  to  105 

0087 

i ndex (k*l) = i ndex (k) 

0088 

105 

k=k  - 1 

0089 

if(t  It  array (k) )  go  to 

100 

0090 

array  (k*l)=t 

0091 

i f (nr  i  1 e  1)  go  to  90 

0092 

i ndex (k ♦ 1 ) sn 

0093 

go  to  90 

0094 

end 
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1 
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0040 
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1 
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0044 
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0045 
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subroutine  wvguid 


This  routine  drives  the 
i nput  el i st . 

If  RPOLY  is  0,  then  a  I  I 
If  RPOLY  is  2,  then  all 
the  routine  RPLYNM. 


generation  of  mode  parameters  using  the 

calculations  are  made  exactly, 
calculations  are  made  approximately  using 


If  RPOLY  is  1,  then  the  initial  calculations  are  approximate  to 
refine  the  initial  solutions  and  the  final  solutions  are  obtained 
using  the  exact  formulation. 


include  'commonl . f or/ I i st ’ 


common/ i nput/f req, rho, az im, cod  r^agf I d , s i gma ,epsr , beta, hpr ime, 

S  hprout 

common/path/path i d , 1 1 ong , 1 1  at , r  '  g , r I  at , rbear , dmax, drm i n, drmax , 

S  year , month, day , gmt, npr i nt, nprof , npath, igcd, ignd,md i r , lost, 

$  I un i t7 ,  I  x 

common/ i onosp/ht I i st (50) , Ini i st (50, 3) , he  I i st (50) , cf I i st (50, 3) , 

S  charge (3) ,mrat i o (3) , nr spec , I htmx , I htmn , I ht , mhtmx , mhtmn ,mht 


character*80  pathid 
integer  year, day 

rea I *4  fr eq,r ho, azim, codip, magf I d , s ■ gma ,epsr , beta , hpr ime, hprout, 

S  1 1 ong, 1 1  at , r I ong, r I  at, rbear , dmax , drmi n , drmax, gmt, 

S  ht I i st , I n I i st, he  I i st, cf I i st,charge,mrat i o 

include  ’common2 . for/I i st ’ 

common /eg  in/el ist(2,30) ,tl ist(2,30) ,dtheta(2) , lub(2) ,deigen(2) , 
S  tht i nc , f  to  I .maxitr, alpha, h,d,prec,»rO,atnmax, debug, typitr, 

S  rpoly.nrtlst 

common/wg  out/tp (30) , tterm(4 , 30) , n term (30) , mode (30) .modes ,  mods 


comple*»8  tp , tterm , dthta 
integer  debug . typ i tr , ^po I y 

rea I »4  e I i st , 1 1 ' st , dtheta , l ub , de i gen , tht i nc , f to  I , a  I pha , h , d , prec , 
S  wrO.atnmax 


equ  .aience  ^ dtheta , dthta ) 
include  common3  for/list’ 

common/ f  f  nc  tn/ omega , • n , theta r , theta i , c , s , esq , ssq , f , df dtht , 

S  hg  , norml 1 , norm?2 , nor  ml ? , rbar 1 1 ,rbar??, 

t  nr  I  ter  ,  neaiq  ,  ad  J  f  I  g  ,  '  sotrp 

rommon/r  matrx/rll,r22,rl?,r?l, 

I  I ogr 1 1  I ogr ?? , I ogr 1? , I ogr ?1 , 

t  dl lldhdi??dh,dll?dh,dl?ldh,ht,delh,topht 

common/m  matrx/ml 1 ,ml?,ml3,m?l , m2? , m23 , m31 , m32 , m33 

' nteger  adj  f I g 

r ea i «8  omega , *r  thetar , theta  <  , ht , de I h , topht , r  f 8 ) , I ogr (8) , d I rdh (8) 
■  >mp ' e « • 1 6  theta  ,  c  s  r sq , ssq , f , df dtht , 
t  hg  norml 1 , nor m2? , norml?, rbar ' 1  r  be  r  ?? 

t  r 11  r??  r 1?, r?l , r s (4) 

1  1 ogr 11,1 ogr  ??  I ogr 1 ? ,  I ogr  ?1  I ogr  * (4)  , 
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0058  1 

0059  1 

0060  1 
0061  1 
0062  1  c 

0063  1 

0064  1 

0065  1 

0066  1  c 

0067  c 

0C68 
0069 
0070 
0071 
0072 
0073 
0074 
0075 

0076  c 

007  7 
0078 
0079 
0080 
0081 
0082 
0083 
0084 
0085 
0086 
0087 
0088 
0089 
0090 
0091 
0092 
0093 
0094 
0095 
0096 
0097 

0098  10 

0099 
0100 

0101  13 

0102 
0103 
0104 
0105 
0106 

0107  15 

0108 
0109 
0110 
0111 
0112 
0113 
0114 


S  dl Ildh,dl22dh,dll2dh,dl21dh,dlrsdh(4) , 

$  mil ,ml2,ml3,m21 ,m22(m23,m31 , m32,m33, 

S  zero/ (0.d0,0.d0)/, zone/ (1 .d0,0.d0)/, 

S  zmp I x i / (0 ■ dO, 1 . dO) / , zdtr / (1 . 745329252d-2 , 0 . dO) / 

equivalence  (thetar , theta) , 

S  (rll , rs) , (I ogrll , logrs) , (d I lldh,d I rsdh) , 

S  (rll,r),  (logrll, logr),  (d 1 1 ldh, d I rdh) 


comp  I  ex* 16  thetaO , stp , rat i o , storel , store2 , store3 , 

$  wterm, ecomp ,m i k 

complex*  8  eigen (30) 
rea I *8  cdang , ref  I ht , capk , stpr , stp i 
integer  psave 
character»20  reason, blank 
equivalence  (el i st, eigen) 

data  blank/’  ’/ , ref  I ht/70 . dO/ 

psave=rpot y 

capk=l  /(I  -  5*alpha*h) 

omega =6  283185306d3*f req 

*n=2  0958426d-2*freq 

wtermzdcmp I x (0  dO, - . 5d0*wn*ref I ht) 

m<  k=dcmp I x (0  dO, -1  d3*»n) 

debug=npr i nt 

adj  f I g*0 

newq*0 

i f (magf Id  le  1  e-10)  then 
i sotrpal 

e  I  se 

if  (cod  ip  eq  90  and  (azim  eq  90  or  azim  eq  270.))  the 
i sotrp=2 
e  I  se 

i SOtrpsO 
end  c  f 
end  i  f 
call  i ntcmp 

>f(rpoly  eq  1)  call  rptynm 
if  (nprmt  gt  0)  print  1010 
kn=0 

ms=0 
i ndex=l 

i f (el i »t(l , i ndex)  eq  0  )  go  to  62 

thetaOse i gen ( i ndex) 

kn*kn*l 

ms=ms* 1 

mnsmode (kn) 

r eason^b I ank 

theta* thetaO 

call  i terat 

fmag*cdabs(f ) 

if(nriter  ge  maxitr  and  fmag  gt  f to  I )  then 
*r i te(reason , 2000)  fmag 
go  to  50 
end  i  f 

pmag*cdab» (r barll*rl2/(zone-rbarll*rll)) 
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0115 

thtr=thetar 

0116 

thti =theta i 

0117 

if(thti  .  ge.  0.)  then 

0118 

wr i te( reason ,2001) 

0119 

go  to  50 

0120 

end  i  f 

0121 

if(kn  .gt.  1)  then 

0122 

do  30  kd=l , kn-1 

0123 

i f (abs (thtr-e 1 i st (1 ,kd))  .gt.  deigen(l)) 

0124 

i  f (abs(thti -el i st(2,kd))  gt.  deigen(2)) 

0125 

write (reason ,2002)  kd 

0126 

go  to  50 

0127 

30 

conti nue 

0128 

end  i  f 

0129 

eigen(kn)=theta 

0130 

33 

i f (ms  .eq.  mn)  go  to  35 

0131 

if(rpoly  .eq.  0  and.  nprint  .gt.  0)  print 

0132 

ms=ms+l 

0133 

go  to  33 

0134 

35 

if(rpoly  .eq.  1)  go  to  60 

0135 

c 

0136 

if  (nr  iter  .gt.  maxitr/2)  then 

0137 

print  *, ’Warning  WVGD:  Excessive  iterat 

0138 

1 ost=2 

0139 

end  i  f 

0140 

s=cds i n  (theta*zdtr) 

0141 

stp=s»capk 

0142 

at=-8 . 6858896d3*wn*d imag(stp) 

0143 

vc— 1  dO/drea 1 (stp) 

0144 

tp(mn)=-zmpl x i *cd log(cdsqrt (zone- stp* stp) ♦; 

0145 

c 

0146 

rat i o=cdsqrt (s) / (df dtht/zdtr) 

0147 

storel=(zone*rbarll)**2* (zone-rbar22»r22) *i 

0148 

store2=(zone»rbarll) * (zone* rba r22) *  rat i o 

0149 

store3= (zone*rbar22) **2* (zone-rbarll*rll) *i 

0150 

ecomp=wterm»storel • (s*hg) **2 

0151 

wm=20  dO«d 1 oglO(cdabs (ecorop) ) 

0152 

wa=cdang (ecomp) 

0153 

if(nprint  gt  0)  print  101 1 , thetaO,mn , nr i 

0154 

S  pmag , at , vc , wm 

0155 

c 

0156 

t  term(l ,mn)=storel/normll 

0157 

t  term(2,mn) =store3/norm22 

0158 

t  term(3,mn)=store2/norml2*r21 

0159 

t  term(4 ,mn)=rl2/r21 

0160 

i f  (cdabs(zone-rll«rbarll)  ge  cdabs(zone- 

0161 

nt eru (mn) =2 

0162 

e  1  se 

0163 

nterm  (mri)  =  1 

0164 

end  i  f 

0165 

go  to  60 

0166 

c 

0167 

50 

if ( r po 1 y  eq  1)  go  to  63 

0168 

if (nprint  gt  0)  print  1012, thetaO, nr i ter 

0169 

if(rho  eq  0  or  npath  eq  1)  then 

0170 

iffrpoly  eq  1)  go  to  63 

0171 

c 

OK  to  drop  a  mode  at  the  transmitter 
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0172 

i f  (kn  .eq.  30  .or.  index  .eq.  30)  then 

0173 

kn=kn-l 

0174 

go  to  62 

0175 

end  if 

0176 

do  53  m=kn,30 

0177 

53 

e  i  gen  (m)  =e  i  gen 

0178 

eiger(30)=(0. ,0.) 

0179 

go  to  13 

0180 

else 

0181 

print  *, ’ERROR  WVGD:  Lost  mode’.mn,’  because  ’.reason 

0182 

!ost=l 

0183 

go  to  999 

0184 

end  i  f 

0185 

c 

0186 

60 

i f  (k  n  .eq.  30  .or.  index  .eq.  30)  go  to  62 

0187 

i ndex=i ndex+1 

0188 

go  to  13 

0189 

62 

nmds=kn 

0190 

eigen(nmds+l)=(0. ,0.) 

0191 

if(nmds  .eq.  0)  go  to  65 

0192 

if(rpoly  .ne.  1)  go  to  999 

0193 

63 

rpoly=0 

0194 

go  to  10 

0195 

65 

print  * , ’ERROR  WVGD:  Lost  all  modes’ 

0196 

lost=l 

0197 

c 

0198 

999 

i f (npr i nt  .gt.  0)  print  1003 

0199 

rpoly=psave 

0200 

return 

0201 

c 

0202 

1003 

format(*  ’) 

0203 

1010 

format(/5x, ’ initial ’,4x, ’mode  i ter ’ ,6x, ’eigen ’ ,8x, ’mag  f’,5x, 

S  'mag  p’ ,5x, ’atten’ ,4x, ’v/c’ ,8x, ’wait’ ’s  exc ’ ,8x, ’theta ’’ ’) 

0204 

0205 

1011 

format(lx,2f7 .3, i 4 , i 5, 2x,2f 7 .3, 2(lx , lpe9 . 3) , lx ,0pf 8 .3, lx , f 9 . 5, 

0206 

S  Ix,f9.3,lx,f6.3,lx,2f8.3) 

0207 

1012 

format (1 x , 2f 7 . 3, 4x , i5,2x,2f7.3,2(lx, lpe9 . 3) ) 

0208 

2000 

format (’fmag= ’ , lpe8 . 2) 

0209 

2001 

format( ’theta i  gt.  O’) 
format(’ it  matches  mode’,i3) 

0210 

2002 

0211 

end 
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